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METHOD AND APPARATUS FOR CREATING A SIMULATED PARTICLE PACK 

BAC KGROUND OF THE INVE NKQNBACKGROUND OF THE INVENTION 
Field of the Invcnti oaField of the Invention 

[0001] The present invention relates to computerized methods and apparatus for 
creating particle packs and, more specifically, to computerized methods and apparatus for 
creating particle packs with enhanced speed and processing efficiency. 

D e scription of th e R e lat e d Art Description of the Related Art 

[0002] The present invention relates to the creation of particle packs that can be used to 
model or simulate particle containing or particle filled systems. There are many varieties of 
particle containing systems. Specific yet merely illustrative examples are found in the fields of 
asphalts, concretes, ceramics, soils, chemical physics of amorphous materials, slurry mechanics, 
and solid propellants or gas generators. An important class of such systems comprises polymeric 
materials with particles embedded or embodied in the polymer body or binder. 

[0003] It is highly desirable in many circumstances to have a capability to model or 
simulate the properties, performance, etc. of such particle containing systems. In designing such 
particle filled systems, for example, substantial trial and error can be avoided if one is able to use 
a model or simulation to assess the implications of making changes in the components, their 
configuration, the conditions, etc. 

[0004] This is particularly true in microstructural analyses, e.g., wherein the analysis 
focuses on the material at a highly magnified level, such as at the molecular or grain stages. A 
specific example of microstructural analysis would involve the modeling or simulation of a 
particle pack for a propellant, which typically would include fuel and oxidizer particles in a 
polymeric binder. When modeling or simulating a particle containing system, it has become 
customary to construct an analytical replica or representation of the particle filled material, 
which is generally referred to as a particle pack. The term "particle pack" has come to be used in 
the field to refer generally to a collection of particles in a defined volume. 
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[0005] The general object in creating particle packs is to model or simulate the manner 
in which a plurality of particles, usually a large number of such particles, will fill or otherwise 
occupy a particular space. The nature of the particles may vary, in some cases considerably, 
depending on the particular application. These variations may include particle size variations, as 
well as other characteristics. 

[0006] The complexity of models used to simulate real world particle containing 
systems increases dramatically if more than a limited number of variables are considered. The 
numbers of interactions, combinations and permutations brought to bear by anything other than 
relatively simple systems causes the processing requirements to become unwieldy even for the 
most advanced computer systems. This complexity and processing requirement limitation is 
even more pronounced for systems that employ large numbers of particles. 

[0007] These limitations have caused workers in the field to invoke various simplifying 
assumptions to render the problem more workable. One such assumption involves the shape of 
the particle. The complexity and related processing demands can be substantially reduced, for 
example, if one assumes that the particles are spherical. 

[0008] The general approach to simulating particle containing systems in the past has 
involved a more rigorous one in which a particle pack is constructed by placing the individual 
particles into a three dimensional volume, i.e., into a "full field." This approach in theory offers 
promise for accurate simulation of the real particle containing systems. The large number of 
particles required to construct a statistically meaningful or useful particle pack, however, usually 
renders this approach impractical. The processing power required to construct such a particle 
pack is prohibitive for all but extremely small numbers of particles, and only relatively small 
variations in particle size. 

[0009] The simplifying assumptions necessary to reduce processing demands down to a 
manageable level can prevent the resulting particle pack from accurately portraying the actual 
packs that are intended to be simulated. 

[0010] An approach that has been used to make the processing demands more 
manageable and yet the accuracy and reliability of the simulation to be good involves a reduced 
dimension approach. An example of this approach is described in I. Lee Davis and Roger G. 
Carter, "Random Particle Packing by Reduced Dimension Algorithms," J.'. Appl. Phy s ., J. 
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Appl. Phys., 67(2), 15 January 1990 pp. 1022-1029 (hereinafter "Davis and Carter 55 ). The 
reduced dimension approach as described there involves simplifying the pack construction 
problem by employing the principle that an arbitrary straight line drawn through a uniformly 
random thr ee- dimensional three-dimensional p article pack (regardless of whether the particles 
are spheres) will, on the average, lie inside particles the same fraction of its length as the packing 
fraction. As the length of the line goes to infinity, the fraction of its length residing inside the 
particles of a random pack goes to the packing fraction exactly. The line is referred to in Davis 
and Carter as a "central string. 55 This approach has proven quite useful in enabling more 
complicated particle packs to be constructed. It can accommodate a substantially wider range of 
particle sizes. 

[0011] Notwithstanding the usefulness of the reduced dimension approach to 
constructing particle packs, its direct application has been limited to a certain extent in that its 
accuracy can be compromised by considering only particles intersected by the central string. In 
real three-dimensional packs, and assuming the particles are spherical, each sphere initially 
intersecting the central sphere, as it slides down the central string during placement, comes to 
rest more often on a "nearest-neighbor 5 ' sphere than on a lower central string particle. The 
sphere then would roll away from the central string, thereby spreading the spheres over a greater 
distance along the central string. 

[0012] A technique for addressing this limitation, also presented in the Davis and 
Carter paper noted above, involves relaxing the requirement that spheres must intersect the 
central. Under this modified approach, particles are permitted to come to rest within a defined 
range of the central string, after assuming that their path is modified to encounter and avoid 
previously placed spheres. This is referred to generally as "perturbation," or perturbed central 
string packing. 

[0013] This modified reduced dimension approach also has been limited, however, in 
that, when perturbation approaches are used, processing demands once again increase, in some 
cases substantially. This can give rise to a fuller and more complicated particle field that once 
again begins to approach the complexity and processing intensity of a full field. 
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Object s of the Invention Obiects of the Invention 

[0014] Accordingly, an object of the present invention is to provide apparatus and 
methods for creating a particle pack that simulate the actual particle distribution within a space 
with a reasonable degree of accuracy. 

[0015] Another object of the invention is to provide apparatus and methods for creating 
a particle pack that involve reduced processing demands and increased processing efficiency 
relatively to full field simulation. 

[0016] Additional objects and advantages of the invention will be set forth in the 
description that follows, and in part will be apparent from the description, or may be learned by 
practice of the invention. The objects and advantages of the invention may be realized and 
obtained by means of the instrumentalities and combinations pointed out in the appended claims. 

SUMMARY OF THE INVENTION SUMMARY OF THE INVENTION 
[0017] To achieve the foregoing objects, and in accordance with the purposes of the 
invention as embodied and broadly described in this document, a machin e impl e m e nt e d 
machine-implemented method is provided for placing a plurality of particles, each having a 
characteristic dimension, to create a particle pack, the plurality of particles comprising N 
categories of the particles, the characteristic dimension of the particles of a given category being 
different from the characteristic dimensions of the particles of other ones of the categories, the 
characteristic dimension of the particles increasing as the category N increases. The method 
compris e s a)d e fining comprises: a) defining a central string, a space disposed about the central 
string, and N concentric subspaces disposed about the central string and within the space, each of 
the N subspaces corresponding to one of the N particle categories; b) selecting a particle from the 
plurality of particles; c) placing the selected particle in the corresponding subspace so that the 
selected particle becomes a placed particle at a particle location unique to that placed particle and 
is in non-overlapping relation with other placed particles, wherein the selected particle plac e ment 
placement, according to one aspect of the inv e ntion invention, includes defining a catch net 
representative of buoyancy of a portion of the placed particles and positioning the catch net 
within the space based upon the placement of the portion of the placed particles. The selected 
particle placement according to another aspect of the invention includes defining a water level 



5 



representative of a level of a portion of the placed particles that are smaller than the selected 
particle and represent a surface of the smaller placed particles, and positioning the water level 
within the space based upon the smaller particle surface. The selected particle-bamg- is placed in 
non-overlapping relation with respect to the catch net and the water level; and d) repeating the 
particle selection (b) and placement (c) until each of the particles of the plurality of particles has 
become one of the placed particles. 

[0018] In accordance with preferred but optional implementations of the invention, the 
following may apply. The particles comprise s phere s spheres and the characteristic dimension 
of each of the particles comprises a radiu s , a radius. The central string comprises-a-fee - a line 
within the space, preferably the z axis of a rectilinear coordinate sy s tem imposed within the 
spacer - z axis of a rectilinear coordinate system imposed within the space. Each of the subspaces 
has a charact e ristic dim e nsion a characteristic dimension that corresponds to the characteristic 
dimension of the particles of the corresponding subspace. Each of the subspaces ma y compris e s 
a cylinder, comprise a cylinder. The subspaces comprise conc e ntric cylinders concentric 
cylinders disposed about the central string, each of the cylinders defining one of the subspaces 
and having a radiu s a radius that corresponds to the radius of the particles in the category 
corresponding to that cylinder. The particle selection may comprise selecting the particles 
randomly or selecting the selected particle from a distribution of th e particles, a distribution of 
the particles. 

[0019] In accordance with another aspect of the invention, an apparatus is provided for 
placing a plurality of particles, each particle having a characteristic dimension, to create a 
particle pack, the plurality of particles comprising N categories of the particles, the characteristic 
dimension of the particles of a given category being different from the characteristic dimensions 
of the particles of other ones of the categories, the characteristic dimension of the particles 
increasing as the category N increases. The apparatus compris e s comprises: a) an input device 
for inputting particle selection information; b) a storage device operatively coupled to the input 
device for storing the particle selection information; c) a processor for selecting a particle from 
the plurality of particles, for placing the selected particle at a particle location unique to the 
selected particle in the corresponding subspace so that the selected particle becomes a placed 
particle at a particle location unique to that placed particle and is in non-overlapping relation 
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with other placed particles, and for establishing a catch net representative of buoyancy of a 
portion of the placed particles and positioning the catch net within the space based upon the 
placement of the portion of the placed particles, the processor placing the selected particle being 
placed in non-overlapping relation with respect to the catch net. 

[0020] In accordance with another aspect of the invention, a machine-readable medium 
is provided for use in placing a plurality of particles, each particle having a characteristic 
dimension, to create a particle pack, the plurality of particles comprising N categories of the 
particles, the characteristic dimension of the particles of a given category being different from 
the characteristic dimensions of the particles of other ones of the categories, the characteristic 
dimension of the particles increasing as the category N increase. The machine-readable medium 
compris e s comprises: a) machine executable instructions for defining a central string, a space 
disposed about the central string, and N concentric subspaces disposed about the central string 
and within the space, each of the N subspaces corresponding to one of the N particle categories; 
b)machin e b) machine executable instructions for selecting a particle from the plurality of 
particles; c) machine executable instructions for placing the selected particle at a particle location 
unique to the selected particle in the corresponding subspace so that the selected particle 
becomes a placed particle at a particle location unique to that placed particle and is in 
non-overlapping relation with other placed particles, the selected particle placement instructions 
including instructions for defining a catch net representative of buoyancy of a portion of the 
placed particles and positioning the catch net within the space based upon the placement of the 
portion of the placed particles, the selected particle being placed in non-overlapping relation with 
respect to the catch net; and d) machine executable instructions for repeating the particle 
selection (b) and placement (c) instructions until each of the particles of the plurality of particles 
has become one of the placed particles. 

[0021] In accordance with another aspect of the invention, an apparatus is provided for 
placing a plurality of particles, each particle having a characteristic dimension, to create a 
particle pack, the plurality of particles comprising N categories of the particles, the characteristic 
dimension of the particles of a given category being different from the characteristic dimensions 
of the particles of other ones of the categories, the characteristic dimension of the particles 
increasing as the category N incr e as e s, th e increases. The apparatus comprises comprises: 
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a) an input device for inputting particle selection information; b) a storage device operatively 
coupled to the input device for storing the particle selection information; c) a processor for 
selecting a particle from the plurality of particles, for placing the selected particle in the 
corresponding subspace so that the selected particle becomes a placed particle at a particle 
location unique to that placed particle and is in non-overlapping relation with other placed 
particles, for establishing a water level representative of a level of a portion of the placed 
particles that are smaller than the selected particle and represent a surface of the smaller placed 
particles, and for positioning the water level within the space based upon the smaller particle 
surface, the selected particle being placed in non-overlapping relation with respect to the water 
level. 

[0022] In accordance with still another aspect of the invention, a machine-readable 
medium is provided for use in placing a plurality of particles, each particle having a 
characteristic dimension, to create a particle pack, the plurality of particles comprising N 
categories of the particles, the characteristic dimension of the particles of a given category being 
different from the characteristic dimensions of the particles of other ones of the categories, the 
characteristic dimension of the particles increasing as the category N increases. The 
machine-readable medium compris e s comprises: a) machine executable instructions for 
defining a central string, a space disposed about the central string, and N concentric subspaces 
disposed about the central string and within the space, each of the N subspaces corresponding to 
one of the N particle categories; b) machine executable instructions for selecting a particle from 
the plurality of particles; c) machine executable instructions for placing the selected particle at a 
particle location unique to the selected particle in the corresponding subspace so that the selected 
particle becomes a placed particle at a particle location unique to that placed particle and is in 
non ov e rlapping non-overlapping relation with other placed particles, the selected particle 
placement instructions including instructions for defining a water level representative of a level 
of a portion of the placed particles that are smaller than the selected particle and represent a 
surface of the smaller placed particles, and positioning the water level within the space based 
upon the smaller particle surface, the selected particle being placed in non-overlapping relation 
with respect to the water level; and d) machine executable instructions for repeating the particle 
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selection (b) and placement (c) instructions until each of the particles of the plurality of particles 
has become one of the placed particles. 



BRIE F DESCRIPTION OF THE DRAWINC SBRIEF DESCRIPTION OF THE 

DRAWINGS 

[0023] The accompanying drawings, which are incorporated in and constitute a part of 
the specification, illustrat e s a illustrate p resently preferred embodiments and methods of the 
invention and, together with the general description given above and the detailed description of 
the preferred embodiments and methods given below, serve to explain the principles of the 
invention. 

[0024] Fig. 1 is a pictorial diagram of a computer system which comprises a presently 
preferred embodiment of the apparatus according to the invention, and upon which the presently 
preferred version of the inventive method may be carried out. 

[0025] Fig. 2 is an illustrative example of a particle pack construction according to the 
preferred implementation, wherein spherical particles are placed along the central string. 

[0026] Fig. 3 shows another illustrative example of a particle pack construction 
according to the preferred implementation, wherein spherical particles are placed along the 
central string, but wherein the particles extend more several particle diameters out from the 
central string. 

[0027] Fig. 4 shows a cylindrical space-and - with cylindrical subspaces within it 
therein according to the preferred embodiment and the preferred method of the invention. 

[0028] Fig. 5 shows an illustrative particle pack constructed using the preferred 
implementation. 

DETAILED DESCRIPTION OF THE DET AILED DESCRIPTION OF THE INVENTION 
PR EFERRED EMBODIMENTS AND METHODS 

[0029] Reference will now be made in detail to the presently preferred embodiments 

and methods of the invention as illustrated in the accompanying drawings, in which like 

reference characters designate like or corresponding parts throughout the drawings. It should be 

noted, however, that the invention in its broader aspects is not limited to the specific details, 

representative devices and methods, and illustrative examples shown and described in this 

9 



section in connection with the preferred embodiments and methods. The invention according to 
its various aspects is particularly pointed out and distinctly claimed in the attached claims read in 
view of this specification, and appropriate equivalents. 

[0030] In accordance with one aspect of the invention, an apparatus is provided for 
creating a particle pack. In accordance with another aspect of the invention, a 
machine-implemented method is provided for placing a plurality of particles to create a particle 
pack. 

[0031] A particle as the term is used herein is used according to its common but broad 
meaning in the field. It includes any type of particle that might be found in a_particle rigid body 
of arbitrary shape containing materials. In the preferred implementation, particles comprise rigid 
or nearly rigid objects of arbitrar y size and shape. 

[0032] It is useful in analyzing particle packs to be able to segregate the particles into 
categories. A group of substantially identical particles, or a group of particles that for purposes 
of analysis may be deemed identical, and which group is in some way distinguishable from other 
groups of particles is called a "mode." If there is only one type of particle comprising the pack 
and those particles are identical to each other, or may be assumed to be identical, the pack is said 
to be "monomodal." If there are two types of particles, the pack is "bimodal," and if multiple 
types of particles are present, it is "multimodal. 55 

[0033] A mode represents a principal size or type of particle. In practice, particles even 
within a mode often are not exactly identical, but instead have slight variations. In some 
instances, for example, each mode can be divided into a size distribution about some mean. The 
distribution may comprise a number of discrete sizes about that mean. The sets of various 
particles of a given size or character within a particular mode are called "submodes." In the 
presently preferred implementation, if there is no size distribution attached to a particular mode, 
the submode arrays simply contain the mode data as their first and only entry. Because both 
modes and submodes theoretically merely involve segregating particles into categories based 
upon one or more distinguishing characters, modes and submodes are referred to generally herein 
as "categories.' 5 

[0034] The apparatus according to the invention preferably comprises a computer 
system 100, such as a commercially available personal computer, small business computer, 
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engineering work station, mainframe, microprocessor, or virtually any othe r machin e s machine 
with similar processing power, appropriately configured and programmed with computer 
software or "code" as outlined and described herein for performing the tasks identified herein, -ft 
The computer system 100 includes a processor 102, storage in the form of RAM 104, a hard 
drive 106, a compact disk (CD) drive 108, and a drive for diskettes 110, all contained within a 
cabinet 1 12. System 100 also includes a monitor 1 16, a keyboard 118, and a pointing device 120 
such as a mouse, track ball or the like. System 100 also may be part of a network, in which case 
it may include a network connection 122 and network bus 124. The method according to the 
present invention, in its presently preferred implementations, may be carried out on computer 
system 100, using the same hardware and software, as described herein with respect to the 
preferred embodiment of the apparatus according to the invention. 

[0035] "Machine" as the term is used herein refers broadly to a computer, processor, 
microprocessor, a network of such devices, or other apparatus capable of performing processing 
as described herein. Such machines typically and preferably will comprise a computer, such as 
computer system 100. 

OVER VIEW OF THE PREFERRED IMPLEMENTA TIO N OVER VIEW OF THE 
PREFERRED EMBODIMENT 

[0036] The method and apparatus according to the presently preferred yet merely 
illustrative versions of the invention will be described in this document as involving the 
implementation and execution of a computer program that can be run on a computer system or 
machine as described above. The preferred versions of the apparatus and methods according to 
the invention are referred to collectively in this document as the "preferred implementation." 

[0037] The architecture of the preferred implementation is outlined below. In its 
presently preferred form, it is written in double precision C++ or Fortran 90. This is not, 
however, limiting. The code may, for example, be written in FORTRAN, or other suitable 
programming languages. An example of such a computer system is shown in Fig. 1. The code 
in its preferred version builds a single particle pack and outputs user-specified pack 
microgeometry and statistics. It is modular and can be disassembled into its component 
subroutines for use elsewhere. 
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[0038] In accordance with the apparatus aspect of the invention, an input device is 
provided for inputting particle selection information, and a storage device operatively coupled to 
the input device for storing the particle selection information. The information used by the 
preferred implementation will be described throughout this document. The input device may 
comprise any one or any combination of several input means, including a network download, the 
CD reader, a diskette reader, the keyboard, the tracking device, a tape drive, or any other means 
used for inputting data into a machine capable of implementing these aspects of the invention. 
The storage device also may comprise any suitable storage means, and may include RAM, the 
hard drive, an optical storage device, storage devices on a machine operably coupled to or 
accessible by computer 100, or other suitable storage devices readable by a computer directly or 
with translation. 

[0039] The main program of the preferred implementation comprises seven subroutine 
calls. The names and functions of these subroutines are as follows: 
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SUBROUTINESUBRO 


SUBROUTINE FUNCTIONSUBROUTINE FUNCTION 


UTINE 

"M A lV/TlT AT A A yf H? 

WAMfcJNAMb 




INITERROR 


Initializes arrays containing the number of errors, warnings, and 
notations accumulated during program execution. It also initializes 
arrays containing information on how these data are to be 
displayed. 


INNAME 


Looks for and reads an input file named INNAME.INP which 
contains the names of files containing the input data needed by the 
preferred implementation to build the user-specified pack and 
output the requested pack information. If the file INNAME.INP 
does not exist, default file names are used to find instructions on 
pack building. Default values are used if user-specified 
instructions are lacking. If essential instructions are lacking, pack 
building is aborted and error messages are issued. 


bb 1 xlllS 1 U(jRAM 


Reads the information for the types, quantities, and sizes of 
particles to be used in constructing the pack. 


SETCONTROL 


Reads flags and parameters controlling the construction of the 
pack. Default values are used if files containing these instructions 
cannot be found. 


SETOUTPUT 


Reads flags and parameters controlling the type of pack data to be 
output. Default values are again used if files containing these 
instructions cannot be found. 


PPACK 


Constructs the particle pack according to instructions read by the 
previous suorouunes. liacn ume u is caneu to duiig a new pacK, it 
initializes all variables and arrays so there will be no interfering 
data left over from a previous call. 


OUTSTAT 


Outputs the user-specified information on pack microgeometry and 
statistics. 



[0040] In summary, the first subroutine (INITERROR) initializes error, warning, and 
notation arrays and parameters. The next four subroutines (INNAME, SETHISTOGRAM, 
SETCONTROL, SETOUTPUT) and the seventh subroutine (OUTSTAT) control the reading of 
building instructions and the output of data. The sixth routine (PPACK) constructs the particle 
pack. This modularity facilitates incorporation of the code or portions of it into larger calling 
codes, such as codes using numerical methods or finite element analysis. 
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[0041] Because in the preferred implementation a substantial amount of information 
may be used to govern particle pack construction, the manner in which the data used by the 
software code is categorized in the context of the preferred implementation will now be 
described. These data categories are placed into common blocks and passed from one code 
section to another as needed. In the presently preferred implementation, all data about the pack 
along with building instructions are placed into one of the following categories: 



DATA BLOCKDATA 


DESCRIPTIONDESCRIPTION 


BLOCK 




HISTO 


Governs the types and quantities of particles to be used in 
constructing the pack, i.e., the particle histogram. 


CONTROL 


Contains all parameters and flags controlling how the pack will be 
DUiii ana sioreu. 


CURINFO 


Contains information on the status of the particle currently being 
placed, such as its position and which other particles it is currently 
touching. 


PLACED 


Contains information on all particles placed so far in the pack 
including the position, size, and type of each particle. 


SETOUT 


Contains the flags and parameters governing what pack statistics 
are to be generated and how the pack microgeometry and its 
statistics are to be output. 


ERRORS 


Contains the names errors, warnings, and notations that can occur 
in program input, execution, and output. It is used to monitor the 
number of times each error, warning, or notation occurs and 
records such information in a file and/or to the screen, according to 
the dictates of the user. 



THE SUBROUTINE INITERROR THE SUBROUTINE INITERROR 

[0042] In the presently preferred implementations, there are three categories of 
information, called "advisors," used to monitor the health of the pack, i.e., to monitor whether 
any of the pack construction rules, such as the prohibitions against particle overlap, and general 
requirements for particle distribution in the pack (e.g., homogeneity at the macro level of the 
random or other specified statistical distributions at the microphone mount level) have been 
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violated. They are contained in the program's error arrays and parameters. The three categories 
of advisors are errors, warnings, and notations. 

[0043] Error advisors are of two kinds, fatal and nonfatal. A fatal error will stop 
program execution because there is not enough information available to build the pack. Nonfatal 
errors will not stop the program execution but will use default values in place of erroneous 
user-data. Therefore, a pack will be built but may not be the kind of pack the user intended. 

[0044] Warning advisors are notifications of unusual situations or potential problems. 
If the warning deals with user input, the user should correct the input problem because the 
warning sometimes precipitates the use of defaults that may not have been intended. 

[0045] The third category is notation advisors. There are many data that may be of 
interest during pack building. Notations report on these items. They are not errors or warnings 
but simply notifications of items of interest as pack building proceeds or after it has finished. 

[0046] All advisors in this implementation have arrays that keep track of how many 
times each has been called; NumErrors(/), NumWarnings(z), and NumNotes(z), for the z-th error, 
warning, or notation, respectively. The totals for each category of advisor are also accumulated 
in the parameters NTotErrors, NTotWarnings, NtotNotes. 

[0047] The display of advisors in this preferred implementation is controlled by arrays 
containing integers limiting the number of times an advisor is to be displayed to the computer 
screen and written to a disk file created for that purpose. The integer arrays LimScErrors(0 and 
LimDiErrors(0 dictate how many times the i-th error is to be written to the screen and the disk 
file, respectively. Similarly, LimScWarnings(z') and LimDiWarnings(/) are provided for 
warnings and LimScNotes(/) and LimDiNotes(z) for notations. If any of these array elements is 
less than or equal to zero, the corresponding advisor is not written. If the user wants all advisors 
written every time they are encountered, the limit arrays need simply be set to a very large 
integer (in this implementation not to exceed 2147483646, which is one less than the maximum 
integer size). The total number of times any advisor has been encountered is monitored in 
NTotErrors, NTotWarnings, and NTotNotes whether or not it is displayed. 

[0048] All data on all advisors according to this preferred implementation are kept in 
and monitored by the subroutine ErrWarnNote(/err, iwarn, inote, iunitnumber), where ierr is the 
integer specifying the error, and similarly for iwarn and inote. The integer iunitnumber is the 
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disk file unit number to which advisors are written as requested. It is assigned the number 13 in 
this implementation. 

THE SUBROUTINE INNAME THE SUBROUTINE INNAME 

[0049] In the presently preferred implementation, the data required by the code are read 
from several data files. The INNAME subroutine reads the names of the data files in which the 
data resides from a file called INNAME.INP. The files whose names are specified in 
INNAME.INP are as follows: 



HLBFILE 


DESCRIPTIONDESCRIPTION 


DEFAULT NAME 






DEFAULT NAME 


Histogram 


Contains data on sizes, quantities, and types of particles 
along with material information and the sizes of the 
cylinders into which they are to be dropped 


HIST.INP 


Drop File 


Contains data on where and how particles are to be 
dropped onto the pack. 


DROPD AT A. INP 


Preplace File 


Containing the positions and types of particles to be 
preplaced in the pack, if any. 


PREPLACE.INP 


Input Control 
File 


Contains the input flags and parameters that control how 
the pack is to be built. 


INFLAG.INP 


Output Control 
File 


Contains flags and parameters controlling what output 
data is to be generated and how it is to be displayed. 


SETOUT.INP 



[0050] In accordance with the preferred implementation, if the file INNAME.INP does 
not exist, then default names are used for each of the five files containing input data. The default 
names for this implementation are identified in the third column of the table above. 

[0051] In the preferred implementation, the histogram file is mandatory but the other 
files are not. Packs can be built if all other files are missing by using defaults. Of course, the 
various options and data production capabilities of the presently preferred implementation may 
not be exercisable without these other files. 
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THE SUBROUTINE SETHISTOCRAM 7?/g SUBROUTINE SETHISTOGRAM 

[0052] The SETHISTOGRAM subroutine according to this preferred implementation 
reads the particle histogram from the histogram file and sets necessary parameters in the program 
associated with that file. The particle histogram comprises a set of particle modes and 
submodes, in this implementation in the form of arrays. If there is no size distribution attached 
to a particular mode, the submode arrays simply contain the mode data as their first and only 
entry. 

[0053] First, subroutine SETHISTOGRAM checks to see if the histogram file exists. If 
it does not, the program is terminated with a fatal error message. If the file exists, the data is 
reads- read as follows. 

[0054] The number on the first line of the histogram file is the "stop build" criterion in 
the form of an integer called STOPBUILD. If STOPBUILD is positive, then that is the total 
number of particles to be placed in the pack. If STOPBUILD is less than or equal to zero, then 
pack building is to be stopped when the center of the highest placed particle reaches a certain 
altitude or when a pr e sp e cifi e d prespecified maximum number of particles is placed, whichever 
comes first. If STOPBUILD is less than or equal to zero, the second line is the real number 
STOP ALTITUDE specifying the altitude at which pack building is to cease. 

[0055] The particle histogram is processed next in this implementation. Each set of 
lines contains the data for one mode. The modes do not have to be in any particular order. The 
first line of the set contains the mode information. If there is a size distribution associated with 
this mode, a set of lines follows that specifies the submodes. In this implementation there are six 
numbers on each mode line. 

[0056] 1 . The first number is the mean particle size for the mode. It is read 
in as a particle diameter but is internally converted to particle radius. 

[0057] 2. The second number is the mass quantity for the mode. Because the 
total mass quantities are normalized later, relative mass quantities are adequate. 

[0058] 3. The third number is the mass density for the mode. The units of 
particle size and density preferably are consistent. In this implementation this is a requirement. 
For example, one should not input a mode diameter in microns and its density in grams/cm 3 . 
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[0059] 4. The fourth number is an integer specifying a material number for 
the mode. Material numbers are a way of specifying some property unique to a mode or set of 
modes. For example, a material number may specify a particular material composition or it may 
specify a set of modes that have some other feature in common. 

[0060] 5. The fifth number is the radius of the cylinder into which the 
particles of this mode will be dropped. The centers of the particles are constrained to stay within 
this cylinder wall. If there is a distribution associated with this mode, the cylinder radius of the 
mode is ignored and the cylinder radii for each submode are used. 

[0061] 6. The last number is a logical parameter RDIST specifying whether 
there is a distribution of submodes associated with this mode. If the parameter is FALSE, there 
is no distribution and the first submode array elements are filled with information for that mode. 
If the parameter is TRUE, the program next reads the submode information in the following 
form. 

[0062] The submode data for this implementation comprises a set of lines, each one of 
which includes the following three pieces of information: 

[0063] 1. The submode size. As with the mode sizes, it is read in as a 
diameter but internally converted to a radius. 

[0064] 2. The fractional mass of the submode. This number is the fraction of 
the mode's mass contained within this submode. The sum of a mode's fractional masses should 
add to unity. The program normalizes them if they do not. 

[0065] 3. The submode cylinder radius. 
[0066] In many parts of the code as presently preferred and implemented, it is 
convenient to have the submodes ordered in increasing particle size. A global size array is 
created that contains all submodes in increasing particle size. The global submode arrays 
comprise an array of submode radii (QSIZE), submode mass (QMASS), submode mass density 
(DENSITY), submode material number (MODEMAT), and submode cylinder radii 
(SUBCYLRAD). 

[0067] The end of a submode is signaled by a set of zeroes everywhere that a nonzero 
number was expected, and the end of a mode is likewise signaled when a set of zeroes is read. 
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[0068] After reading all mode and submode information, the SETHISTOGRAM 
subroutine reads the information in the drop file if it exists. If it does not exist, the program 
drops the particles randomly anywhere within the cylinder corresponding to that mode or 
submode. If it does exist, then it contains the following data. 

[0069] The first line has two integers. The first integer is a drop option, either 1, 2, 
or 3. The second integer is the number of particles that are to be dropped according to that 
option. 

[0070] Option 1 means that first group of particles are to be dropped randomly 
anywhere within their cylinder. 

[0071] Option 2 means that the particle group is to be dropped within a circle of 
relative radius a r where the drop radius is some fraction of the cylinder radius (0 < a r < 1) of the 
particle that will be selected. The center of this circular region is at relative radial position p r 
where p r is some fraction of the particle's cylinder radius when it is chosen (0 <p r < 1). Because 
dropping particles outside the cylinder wall is forbidden, it should hold that p r + a r < 1 . These 
data are stored in the pack position arrays. 

[0072] Option 3 means that each particle in that group is to be dropped from a specified 
position and is of a specified mode. The x position of the drop site is stored in the x position 
array, the y position in the y position array, and the mode number in the mode number pack 
array. If an-fe^V (x % v) p osition lies outside the cylinder corresponding to that mode, a warning 
message is issued and the particle is randomly dropped anywhere within the cylinder. 

[0073] The end of this input file is signaled by a value of zero for the drop option. 

[0074] Another function of the SETHISTOGRAM subroutine in this implementation is 
to create an index array that orders the global cylinder radii in ascending order with an integer 
array called CylSizeOrder so that cylrad(CylSizeOrder(l)) is the smallest cylinder radius, 
cylrad(CylSizeOrder(2)) is the second smallest, and so forth. 

THE SUBROUTINE SETCONTROL THE SUBROUTINE SETCONTROL 

[0075] The currently preferred implementation includes a variety of options for pack 
construction. Some of those options are read in a file that can have any name as specified in 
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Subroutine INNAME but whose default value is INFLAG.INP. This subroutine SETCONTROL 
reads those control parameters. 

[0076] If a data file containing input control parameters exists, then the control 
parameters are read from it in this implementation. If it does not exist, a default is set for each 
parameter in this routine. The parameters used in this implementation are described in detail 
below but in summary include the following: 

[0077] ROUGHNET = A flag that specifies whether the^-fe - k-th mode catch net 
(described below) is ratcheted up after each particle placement so that no two particles are at 
exactly the same altitude. It is mostly used at the beginning of pack building so the spheres at 
the bottom of the pack do not lie exactly on the bottom of the cylinder. The default is TRUE. 

[0078] POOL = A flag that specifies whether the number of particles in the pack are 
selected by a random number generator, weighted by their frequency of occurrence in an infinite 
pack, or whether they are selected from a pool whose mode proportions are, to the nearest 
integer, exactly the proportions of an infinite pack. If the POOL flag is TRUE, then random 
chance variations in modes that contribute few particles to the pack are avoided. The default is 
TRUE. 

[0079] FINDLOW = A flag that orders find the lowest exposed north pole of placed 
particles within its drop cylinder and to drop immediately in its neighborhood. It is overridden 
by any instructions that may appear in the drop file. The default is TRUE. 

[0080] NDROP = A positive integer that specifies the number of times a trial drop of 
each particle is to be made. The lowest final resting place is selected from among the trial drops. 
Higher values for NDROP means the pack will be more dense. The default is one. 

[0081] ISEED = An integer specifying the random number seed. A non-positive 
integer will trigger to seed the random number generator from the system clock. 

[0082] TUNNEL = A flag that gives permission to check for cavities underneath larger 
spheres when a smaller sphere is dropped. The small sphere "tunnels" through the large one and, 
if there are no overlaps in the southern or underneath side, is dropped again from that underneath 
position whereupon it searches for a final resting place in that cavity. The tunneling procedure 
seeks to mimic tamping or shaking a pack to fill in such cavities, thereby making the pack more 
dense. 
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THE SUBROUTINE SETOUTPU T THE SUBROUTINE SETOUTPUT 

[0083] A particle pack is typically created so that the system or method user can obtain 
specific information about the actual particle pack the simulation is intended to simulate. The 
information required depends on the application. The subroutine SETOUTPUT in this 
implementation reads the flags and parameters that instruct other parts of the preferred 
implementation to create specific pack information. In this implementation, for example, if the 
instruction file for setting output flags and parameters does not exist, the code outputs a limited 
amount of information into a file called PACKDATA.OUT. Otherwise, this basic pack 
information is output to a file name specified in Subroutine INNAME. With no specifications, 
the following data is sent to the output file: 

• an echo of the particle histogram; 

• the packing fraction; 

• the maximum and minimum positions of particles for each submode as 
well as overall maxima and minima; and 

• the number of particles placed for each submode as well as total particles 
placed. 

[0084] If the pack output file exists when using the preferred implementation, it 
preferably contains the following information. The first line contains a flag specifying whether 
the particle pack is to be written to a disk file. The default is FALSE. If the flag is TRUE, then 
the next set of lines sp e cify specifies how the pack file is to be written. On the first subsequent 
line, the format of the pack (ASCII or binary) is written. On the second subsequent line, a flag is 
read that determines whether the material number and radius are written in the pack file in 
addition to the three coordinates and the mode number. On the third subsequent line, the output 
geometry (Cartesian or cylindrical coordinates) is written. The name of the file to which the 
pack is to be written is on the fourth line. 

[0085] The next line contains a flag specifying whether radial distribution functions 
gi/r) are to be generated. If it is TRUE, then the functions that are to be generated are specified 
in a list of integer pairs i and j on the subsequent lines, where the / and j represent modes or 
submodes. Each line also contains a Ar which is the increment in terms of the y'-th particle radius 
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in which the number of particle centers are to be counted (default = 0.1 radius) and r max which is 
the extent to which the radial distribution function extends in terms of y'-th particle radius. Hence 
each radial distribution function requested is on a line with the format 

* j Ar r max . 
The list ends when zeroes are read. Each radial distribution function is written to a file with the 
name g[i]_D']-dat where-the - each set of square brackets indicat e the indicates an actual specified 
integ e rs integer to be used in the file name, e.g., g3_l 4.dat for the radial distribution function of 
Mode 14 spheres about Mode 3 spheres. The default for the radial distribution flag in this 
implementation is FALSE. 

[0086] The next line in this implementation contains a flag specifying whether 
coordination numbers are to be generated. If it is TRUE, then the coordination numbers that are 
to be generated are read as a list of integer pairs / and j on subsequent lines. The list ends when 
zeroes are read. The coordination numbers preferably are all written to one disk file 
COORNUM.DAT. The default for the coordination number flag in this implementation is 
FALSE. 

[0087] Radial distribution functions and coordination numbers and their calculation 
procedures are defined in the discussion of subroutine OUTSTAT, below. 

[0088] In accordance with the apparatus aspect of the invention, the apparatus 
comprises means for defining a space, a plurality of subspaces within the space, and a central 
string within the space and within the subspace. In accordance with the method aspect of the 
invention, the method comprises defining a space, defining a plurality of subspaces within the 
space, and defining a central string within the space and within the subspaces. 

[0089] The subspaces preferably comprise N concentric cylinders disposed about the 
central string and within the space. Each of the N cylinders corresponds to one of the N particle 
modes and has a cylinder radius corresponding to the particle radius of the particles of the 
corresponding particle mode. The means as implemented in the preferred embodiment 
comprises the processor. 

[0090] The space and subspaces are constructions that define the regions or volumes 
into which the respective particles are to be placed. The space may comprise any contained 
geometric region, but preferably comprises a cylindrical construction. The subspaces also may 
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comprise any contained geometric region, but preferably have a shape and configuration that 
corresponds to the overall space. In the preferred versions of the apparatus and method, the 
space comprises a cylindrical volume, and the subspaces comprise a plurality of concentric 
cylinders, disposed about a longitudinal or z axis, wherein sequential ones of the cylinders have a 
radius that is larger than the prior cylinders. An illustrative example of such a concentric 
cylinder configuration is shown in Fig. 2. Fig. 2. The number of cylinders in the preferred 
implementation is equal to the number of particle categories, e.g., submodes, that will be present 
in the particle pack. The radius of each cylinder is selected to correspond to the radii of the 
particles. The radius of the first cylinder corresponds to or is a function of the particle radius of 
the smallest mode, the radius of the second cylinder is equal to the particle radius of the second 
smallest mode, and so on. The cylinder radius for a particular particle category preferably 
corresponds to or is a function of the radius or a representative radius of the corresponding 
particle for that mode or submode. The cylinder radius may be any positive number multiple of 
the particle radius. Pr e f e rably but optionally Preferably, but optionally, the cylinder radius is 
about 3 to 10 times the particle radius, and mov e pr e f e rably more preferably, about 5 to 10 times 
the particle radius. 

THE SUBROUTINE PPAC K THE SUBROUTINE PPACK 

[0091] Further in accordance with the method aspect of the invention, the method 
includes selecting a particle from the plurality of particles, and placing the selected particle in the 
corresponding subspace so that the selected particle becomes a placed particle at a particle 
location unique to that placed particle and is in non-overlapping relation with other placed 
particles. 

[0092] In accordance with the preferred implementation, particles are selected 
randomly from the distribution of particles using a Monte Carlo method, although this is not 
necessarily limiting. In terms of the apparatus, this comprises a random number generator in 
computer system 100. 

[0093] Regarding terminology, a "selected particle," also referred to herein as a 
"current particle," is a particle that has been selected for placement, but has not yet been placed 
(become a plac e d "placed p article") in a particle location in the particle pack. A "placed 
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particle" is a particle that has been placed or fixed by assigning it a unique particle location in the 
particle pack. 

[0094] In the preferred implementation, the software code subroutine selects particles 
and constructs the particle pack according to instructions read by the previous subroutines. Each 
time it is called to build a new pack, it initializes all variables and arrays so there will be no 
interfering data left over from a previous call. Code lines of PPACK and a brief explanation of 
each will now be provided. 

CALL PREPAREPAC K CALL PREPAREPACK 

[0095] Subroutine PREPAREPACK performs three functions in this implementation. 
First, it preplaces particles, if any, that have a position assigned by the user, such as would be 
desirable i f s e mi random semirandom p acks were to be built. Secondly, it arranges histogram 
data such that particles can be chosen by the random number generator with probabilities 
appropriate for their abundance in the reduced dimension pack. Third, it initializes and sets 
parameters that will be used by the rest of the subroutines in PPACK. 

DO icarus = 1, stopbuild 

[0096] This DO loop continues to cycle until either the user-specified number of 
particles-have - has been placed in the pack or until a pack height has been reached. 

CALL CHOOSE CALL CHOOSE 

[0097] Subroutine CHOOSE chooses the type, category, mode or submode of the next 
particle to be dropped. It drops the "current particle" according to user sp e cifi e d user-specified 
instructions. 

DO idrop = 1, ndrop 

[0098] The user may specify more than one drop at different drop sites for each particle 
so that the lowest final site can be selected. Multiple drops are equivalent to tamping the pack so 
the particles settle into lower sites resulting in a denser pack. 
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CALL DROPSIT E CALL DROPSITE 

[0099] Subroutine DROPSITE selects an-tf^)- (% Y) p osition from which the particle 
will be dropped. The site can be chosen randomly, or the user may ask the code to use a criterion 
to find potential sites where the particle is likely to settle to an especially low final resting place. 

CALL PLACE CALL PLACE 

[00100] Subroutine PLACE finds the final resting place or "placed particle location" for 
the current particle, which is the site at which the current particle will ultimately be placed. The 
subroutine PLACE in this implementation contains all the instructions on how a particle drops, 
rolls, and encounters other objects in the pack as it moves toward its final resting place. It is 
where the machine spends the great majority of its processing time during pack construction. 

END DO 

CALL UPDAT E CALL UPDATE 

[00101] Subroutine UPDATE updates all the particle pack arrays and parameters with 
the position and additional data associated with the newly placed particle. 

END DO 
RETURN 
END 

[00102] In accordance with the preferred implementation, a particle pack is constructed 
by "dropping" the selected particles one at a time into the set of concentric cylinders. The 
dropping presumes a vertical orientation, which is not necessarily required. But regardless of the 
reference frame, the current particle is moved in the space until it encounters another object, as 
described in principle herein. Each particle mode is dropped somewhere within the cylinder 
radius of its corresponding cylinder. The radius of the cylinders in this preferred implementation 
is specified by the user, but preferably the small particles are dropped into a small-radius 
cylinder and the larger particles into a larger-radius cylinder or subspace. The large particles can 
therefore overlap the inne r small radius small-radius cylinders but the small particles do not 
have to fill the volume in the larger cylinders. The reduced dimension approach as described in 
Davis and Carter is used in this preferred implementation. 
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[00103] The reduced dimensional approach can make pack building tractable from a 
processing demand perspective when the particles comprising the pack have broad size ranges. 
Particle centers in this implementation cannot pass outside the radius of the cylinder wall to 
which they belong but can pass through any smaller cylinder walls. Because only the innermost 
cylinder can have all modes or submodes intersecting it, only there does the reduced dimension 
pack literally approximate the corresponding real three-dimensional pack. The remainder of the 
pack can however be simulated using this approach, for example, by assuming that the remainder 
of the space is statistically identically configured. 

[00104] In the presently preferred implementation, the selection of particles to be placed 
and the placement of the particles in the particle pack are undertaken by the processor within the 
machine under the execution of the PPACK subroutine. In this implementation, there are four 
subroutines called by PPACK (PREPAREPACK, CHOOSE, PLACE, and UPDATE). Each of 
these subroutines and the functions they perform will now be described. 

THE SUBROUTINE PREPAREPAC K THE SUBROUTINE PREPAREPACK 

[00105] As noted above, the PREPAREPACK subroutine has three functions: 
preplacing particles, setting cumulative arrays so the random number generator can readily 
choose the mode of the next particle to be placed, and initializing various parameters in 
preparation for building the pack. Each of these functions will now be described in turn. 

Pre Placing Particl e s Preplacing Particles 

[00106] The first function of PREPAREPACK is to pr e plac e preplace any particles for 
which the user has specified a position. This option permits s e mi random semirandom or fully 
ordered packs to be constructed. A filename is read in Subroutine INNAME that contains the 
global (X, 7, Z) position and the mode or submode number of each particle to be pr e plac e d. 
preplaced. If the filename is not read, the default filename for preplaced particles is 
PREPLACE.INP. If a filename of preplaced particles does not exist, the program assumes no 
particles are to be preplaced and this option is not exercised. The preplace file contains five 
numbers for each particle to be placed. In order, they are the X position, the Y position, the Z 
position, and the mode number, and the submode number. These are read into the surface arrays. 
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Setting Arrays For Choo s ing Particlcs Setting Arrays for Choosing Particles 
[00107] The second function of PREPAREPACK according to the preferred 
implementation is to make an array so that the random number generator can randomly choose, 
with the appropriate probability, the submode of the next particle to be placed in or on the pack. 
The relative fraction of particle numbers for each submode therefore should be calculated. The 
number of particles in each submode to be placed will depend on the submode mass, its mass 
density, its particle radius, and its cylinder radius, e.g., in the following manner. 

[00108] In this implementation there are two options in setting up the weighting factors 
that govern how the random number generator will choose the next particle to be placed. The 
first option is to use the same weighting factors each time a particle is chosen for placement. 
This is referred to herein as the "unchang e d w e ighting "unchanged-weighting option." The 
second option, referred to herein as a "pool option," creates a pool that has a set of bins with the 
appropriate number of particles placed into it. Each time a particle is drawn from one of the 
bins, the weighting factor of drawing from each bin is altered slightly. Creating a pool of 
particles from which to draw forces the right proportion (to within nearest integers) of particles 
from each submode to be placed. In the limit of an infinite number of particles in the pack, these 
two methods are identical. The pool option more closely resembles a very large pack. The 
unchanged-weighting option will be described first. 



Unchanged - Weighting Pool Oprifl« U nchanged- Weighting Pool Option 
[00109] The total volume of A>th submode particles to be placed is 



m. 



4 3 



(1) 



where m k is the submode mass, p k is the submode mass density, and a k is the particle radius for 
that submode. 

[00110] The relative number of submode k particles in this reduced dimension pack is 
also proportional to the area of the confining cylinder for that submode: nR k where R k is the 
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radius of the cylinder. Therefore, the fraction fk of £-th submode spheres to be placed in a 
reduced dimension pack is proportional to 



<k = 



m k R 2 k 
Pk<*l 



(2) 



[00111] To obtain the actual fraction fk, the above expression can be normalized: 



(3) 



where TV is the total number of particle submodes. Then 



fk = 9JF. 



(4) 



These are the unchanged-weighting factors that determine which particle submode shall be 
chosen for the next placement. 

[00112] These weighting factors are placed into a cumulative array so the random 
number generator can immediately choose the submode. In other words, the cumulative array is 
the sum of all weighting factors prior to it: 



[00113] When a uniform random number r is generated on |"0J"|, [0, 1], a bisection 
search is made of the arrav-f^T^TTrrGAri - iCu C?, . . . Gyl to find which pair of numbers the 
random number lies between. The k-th submode is chosen if Q.y < r < C*. Co is defined to be 
zero. 



C f =£ ft, \ <i<N 



(5) 



C N = \. 
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Pool Option PooX Option 

[00114] The pool option is similar. As discussed before, the user can specify the 
number of particles to be placed in the pack or can specify a pack height in which case the 
program continues to add particles until the height is achieved (or until the array sizes are 
exhausted). The pool option in this implementation can only be used if the user has specified the 
total number of particles to be placed in the pack, say n. The number of submode k particles nk 
to be placed is 

n k = NINT[«/*] (6) 

where NINT indicates the integer nearest to its argument. Because of nearest integer roundoff, it 
is possible that the rtk do not sum exactly to n. Then the rtk quantities should be adjusted to 
compensate where the roundoff errors were the largest. 

[00115] It is possible, especially for packs with very broad particle size distributions, 
that some nf k < 14 so that submode would not be represented at all in the pack. The preferred 
implementation issues a warning when that happens, or even when a particular n k is very small, 
so that the user is aware that the pack is too small to be able to sample enough different 
configurations to mimic a true three-dimensional particle pack. Changing cylinder sizes on some 
submodes can solve this situation, should it arise. 

INITIALIZING VARIOUS PAILiMETER S hiitmlizmz Various Parameters 

[00116] The packing routine can be called many times, as when the user wishes to 

search for an optimum pack subject to some set of constraints. Some of the parameters and 

arrays from previous calls must not be carried into successive calls or they can taint the results. 

This portion of Subroutine PREPAREPACK initializes all necessary parameters, flags, and 

arrays for use in a new pack simulation. 
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THE SUBROUTINE CHOOS E THE SUBROUTINE CHOOSE 

[00117] The CHOOSE subroutine uses the random number generator to choose the 
submode of the next particle to be added to the pack. The random number generator in this 
implementation generates a uniform distribution on [0J1> [0, 1]. The cumulative array from the 
previous routine PREPAREPACK is a set of numbers such that 



[00118] Because this cumulative array can be large, and because there is no restriction 
on the relative distances between successive array elements, in this implementation, a search by 
bisection is a preferred method to find two numbers between which a given random number lies. 

[00119] When a random number r is chosen in \0A], |U 1] we want to find the two 
array elements that bracket it such that 



(In the unlikely event that r = 0, it will be assumed to belong to the first bin, i.e., where k = 1.) 
In preferred bisection methods, a check is made to determine on which side of Cm/ 2 the random 
number lies. Then attention is concentrated on that side and it is bisected again, continually 
narrowing the range of {Ck} in which r lies until Eq.-(8)-_8_is satisfied. 

[00120] As noted above, the particle just chosen is referred to herein as the "current 
particle" "current particle" or the "selected particle" in the discussion below. 

THE SUBROUTINE DROPSIT E THE SUBROUTINE DROPSITE 

[00121] The DROPSITE subroutine chooses a drop site from which the current particle 
can be dropped. Actually, if the user so specifies, it chooses several drop sites and then keeps 
the site at which the particle came to the lowest possible final resting place, referred to herein as 
a "placed particle location." In a sense, this is the same as tamping the pack so that it settles into 
a denser configuration. The user specifies how many trial drops a particle is to undergo. The 
default in this implementation is one. 



Co = 0<Ci<C2<... C N .j<C N =1. 



(7) 



C k ,j<r<C k ,l<k<N 



(8) 
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[00122] There are several options available to the user in this implementation when 
deciding where a new selected particle is to be dropped onto the pack. 

[00123] 1 . Drop it randomly anywhere within its cylinder; 

[00124] 2. Drop it randomly anywhere within a user-specified circle that is 
contained within the cylinder; 

[00125] 3. Drop it from a user-specified-fe^j^ - (X* Y) p osition; and 

[00126] 4. Select the lowest exposed north pole among the placed pack 
spheres and drop it randomly within a short distance of that pack sphere' s-fe¥^) - (X, Y) p osition. 
If dropping multiple times, it is preferable to select the several lowest exposed north poles and 
drop within a short distance of each of them. 

Dropping Randomly Within The Cv/t/itfe r Dropping Randomly Within the Cylinder 
[00127] When dropping a sphere anywhere within a given cylinder, two random 

numbers are generated, u 9 and u p . The azimuthal angle (in radians) of the drop position <p c 

relative to the central string is chosen by 



[00128] The radial position is not chosen uniformly because there is more area for large 
radial positions in a circle than for small radial position. Because the area grows as the square of 
the radial size, all areas will be equally accessible if we take the square root of the random 
number generator and multiply it by the cylinder wall radius for the current sphere W c \ 



(9) 



Pc = W c Ju~ p . 



(10) 



[00129] The (X, Y) drop site is found as follows: 



X c = p c cos <p c 



(11a) 



Y c = p c sin <p, 



c • 



(lib) 
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Dropping Randomly Within A User - Specified CfVcte Dropping Randomly Within a 
User-Specified Circle 

[00130] Suppose the user specifies that the particles are to be dropped randomly within a 
circle of radius R centered at (Xq 9 Yo). In the preferred implementation it is required that the 
circle be contained within the cylinder of the current sphere, i.e., 

y/Xl +yg +R<W C . (12) 

Otherwise, the circle is an illegal specification because a particle cannot be dropped outside its 
cylinder wall. In that case, the preferred implementation will default to dropping randomly 
anywhere within the cylinder wall at W c . A value of zero for R is permitted, meaning the 
particles are dropped from the position (Xo, Yo). If this condition is satisfied, then the drop site 
can be selected as follows. Two random numbers are generated, u<p and w p . Then 

(pc^lnuy (13) 

Pc = R^ P (14) 
X c =X 0 + p c cos <p c (15a) 



and 



Y c = Y 0 +p c sin q> c . (15b) 

Dropping From A User - Specified Xy Position Droppim From a User-Specified XY 
Position 

[00131] The condition of dropping a selected particle from a user-specified X-Y position 
is a subset of the previous discussion for the special case of R = 0. The positions are stored in the 
pack surface arrays for that particle until it is placed, at which time the drop XY position is 
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replaced by the XY position of the final resting location, or the placed particle location, for that 
particle. 

Dropping Above Lowest iS/fcs Dropping Above Lowest Sites 

[00132] In the preferred implementation, a running list is kept of all placed pack spheres 
whose north poles have not yet been covered by other pack spheres. These north poles are 
ordered according to their altitude (north pole altitude, not sphere center altitude). It is generally 
desirable for the current or selected sphere to find a final resting place with the lowest altitude. 
The lowest exposed north pole is a more likely candidate for a low final placed particle location. 
Therefore, if the current sphere is to be dropped m times (m > 1), then we choose the lowest m 
exposed north poles and drop randomly within a circle of radius MIN(ot c , a t ) of that north pole's 
XY position, as long as it is within the current sphere's corresponding cylinder wall. The 
quantity a c is the radius of the current sphere and a, is the radius of the placed sphere above 
whose exposed north pole the current sphere is being dropped. Dropping randomly within this 
circle is done in the same manner as is described above. North poles outside W c are ineligible for 
dropping in this implementation. If there are less than m exposed north poles, the other sites are 
chosen randomly within the cylinder. 

THE SUBROUTINE PLACE THE SUBROUTINE PLACE 

[00133] The subroutine PLACE in this implementation takes the current particle from its 
initial drop site to the final resting place or placed particle location associated with that drop site. 
The dropping and rolling preferably is done without any dynamics. The particles in this 
implementation do not bounce, nor do they gain speed or shoot off the edge of anothe r particl e s 
particle with a parabolic trajectory. It is as if the particles are being dropped in a highly viscous 
fluid so that all inertia is absent and the particle creeps to its final resting place. 

[00134] There are four subroutines in the presently preferred implementation of the 
PLACE subroutine. They are referred to herein as FOBJ1, FOBJ2, FOBJ3 and STABLE. 
FOBJ1 finds the first object the descending current sphere encounters. FOJB2 finds the second. 
FOJB3 finds the third. STABLE determines which contacted objects the current sphere will stay 
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in contact with and which it will roll away from. The term "object" is used instead of "sphere" 
because one of the objects of contact may be the cylinder base or wall. 

[00135] A number of terms will now be defined for purposes of this discussion before 
describing the operation of the PLACE subroutine. The "roll corridor" is the path of descent of a 
sphere as it descends along its contacts with one or more spheres and/or the cylinder wall. A 
good share of the effort of subroutine PLACE is to find the objects that overlap the roll corridor. 
An "object" according to this implementation can be one of four things that the current sphere 
can encounter as it descends: (1) other pack spheres (placed particles), (2) the cylinder base or 
wall, (3) a "water level," or (4) a "catch net." The latter two are described more fully below. 
"Contact stability" refers to whether the current sphere is in compressive contact with an object 
or in tensile contact. If it is in compressive contact with another sphere or the cylinder wall, they 
are pushing against each other due to the weight of the current sphere when gravity is acting 
along the -z axis. The current sphere at this position stays in contact with an object it 
compressively touches. Usually if it touches three or more objects compressively, it is stable and 
is placed at that position. If there is a tensile contact, the object and current sphere are pulling 
away from each other. Therefore, the current sphere will roll away from the object unless it 
already has three or more compressive contacts. The first pack sphere touched by the current 
sphere is denoted as Sphere 1, the second as Sphere 2, and so forth. If a set or plurality of 
spheres-are- is touched simultaneously, the numerical order can be arbitrarily selected. 

[00136] Th e flowchart pr e s e nted in Fig. shows how Subroutine PLACE op e rat e s. 
Subroutine FOBJ1 searches for the first object of contact. If the current sphere simultaneously 
contacts more than one object, it checks them for stability through subroutine STABLE. 
Subroutine FOBJ2 searches for the second object of contact. If the current sphere 
simultaneously contacts two or more objects, it checks them for stability. Likewise, subroutine 
FOJB3 searches for the third object of contact. Again, if the current sphere simultaneously 
contacts three or more objects, it checks them for stability. The current sphere is placed if it is 
stably touching three or more objects or has hit the catch net or a water level. 

[00137] We will now give an overview of each subroutine and how they interact, 
followed by a detailed discussion of each one. 
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Overview of Subroutine STABLE Overviev/ of Subroutine STABLE 
[00138] Subroutine STABLE determines whether the current sphere is in compressive or 
tensile contact with all objects it is currently touching. It marks those objects in tensile contact 
so that the current sphere will roll away from them as it continues to descend. 

Overview of Subroutine FOBJ1 Overview of Subroutine FOBJ1 

[00139] The first subroutine for determining object contact is called FOBJ1 for "Find 
Object 1," object 1 being which is the first object the descending current sphere will encount e r. 
encounters. There are four possible outcomes o f FOBJ1; FOBJ1: 

[00140] • 1. The current sphere encounters the catch net or a water level 
before encountering any other sphere. If it does, then it is placed at the position where such 
contact occurs. 

[00141] • 2. The current sphere encounters a placed sphere of the pack. 
If it does, the index number of that sphere (hereafter denoted as Sphere 1) is recorded and FOBJ2 
is called. 

[00142] • 3. The current sphere simultaneously hits two placed spheres 
of the pack. If it compressively touches both of them, it will roll along the roll corridor defined 
by them and FOBJ3 is called to determine whether a third object is hit before it loses contact 
with one or both of the touched placed spheres. 

[00143] • 4. The current sphere simultaneously hits three or more placed 
spheres of the pack. STABLE checks to see if any three of them provide a stable placement. If 
they do, the current sphere is placed there. If they do not, STABLE marks the objects from 
which the current sphere will roll away. 

Overview of Subroutine FOA/i Overview of Subroutine FOB J2 

[00144] FOBJ2 searches for the second object to be encountered by the current sphere as 
it descends along the path of steepest descent on Sphere 1 . There are four possible outcomes of 
FOJB2: 

[00145] • 1 . The current sphere hits the catch net or a water level before 
encountering any other object. If it does, it is placed there. 
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[00146] • 2. The current sphere loses contact with Sphere 1 by reaching 
its equator before encountering any other object. If it does, it is free falling again and FOBJ1 is 
called, and re-initialized. 

[00147] • 3. The current sphere hits a second object, either another 
placed sphere or its cylinder wall. STABLE is then called to see if it is in compressive contact 
with both objects. If it is, FOBJ3 is called. If it loses contact with one object, FOBJ2 is r e call e d 
recalled and re-initialized. If it loses contact with both objects, it is in free fall again and FOBJ1 
is r e call e d recalled and re-initialized. 

[00148] • 4. The current sphere hits two or more additional objects 
simultaneously. STABLE is called to check for compressive contacts. If there are three or more 
compressive contacts, the particle is placed where it is. If there are two, FOBJ3 is r e call e d; 
recalled; if there is one, FOBJ2 is r e call e d; recalled; and if there are none, the current sphere is 
in free fall and FOBJ1 is re called, recalled. 

Overview of Subroutine FOBJ3 Qverview of Subroutine FOB J3 

[00149] FOBJ3 searches for the third object the current sphere encounters as it descends 
along the roll corridor formed by the first two objects. There are four possible outcomes of 
FOBJ3. FOBJ3: 

[00150] • 1. The current sphere hits the catch net or a water level before 
encountering any other object. If it does, it is placed there. 

[00151] • 2. The current sphere loses contact with one of its two 
contacts before encountering any other object. If it does, FOBJ2 is r e called, recalled. 

[00152] • 3. The current sphere loses contact with both objects 
simultaneously before encountering any other object. If it does, it is in free fall and FOBJ1 is 
r e call e d, recalled. 

[00153] • 4. The current sphere encounters one or more additional 
objects. If it does, STABLE is called to check whether three or more are in compressive contact. 
If they are, the current sphere is placed there. If not, contacts are lost with those objects in 
tensile contact. If two objects remain in tensile contact, FOBJ3 is r e called; recalled; if one 
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object, then FOBJ2 is r e call e d; recalled; and if all contact is lost, the current sphere is in free 
fall and FOBJ1 is r e call e d recalled again. 

Detail of Subroutine STABL E DeizxX of Subroutine STABLE 

[00154] The structure and manner of operation of the STABLE subroutine according to 
the preferred implementation will now be addressed in greater detail. Consid e r Particularly, the 
STABLE routine may consider the current sphere with weight W acting in the -z dir e ction. If 
direction, if it is touching one or more obiects,4s4t - if it is compressive! y touching th e m or would 
it - them, or if it would roll away from-them r the one or more objects. 

[00155] Index numbers of placed pack spheres being touched are passed to subroutine 
STABLE in an array called INDEXTOUCH. A separate parameter called TOUCHWALL states 
whether the cylinder wall of the current sphere's submode is being touched. The manner in 
which the preferred implementation determines the compressive touching criteria when the 
current sphere touches (a) one sphere, (b) two spheres, (c) three spheres, (d) one sphere and the 
wall, or (e) two spheres and the wall will now be addressed. There is generally no need to 
consider compressive touching for more than three objects because the case for more than three 
objects can be broken into sets of three. At least three objects should be checked, however, 
because that is the minimum number needed for stable placement in the presently preferred 
implementation. 

Compressive touching of one sphere Compressive touching of one sphere 
[00156] Determining compressive touching on one and only one pack sphere can be 
done as follows. If the current sphere is touching its northern hemisphere, it is compressive. 
Because current spheres are dropped from above, they cannot initially touch southern 
hemispheres and FOBJ1 rejects a sphere as a touching object if the current sphere only touches 
its equator. 

Compressive touching of two spheres Compressive touching of two spheres 
[00157] Consider a current sphere that has descended along a longitud e longitudinal 
line of Sphere 1 until it reaches Sphere 2. When it touches Sphere 2, there are only two 
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possibilitios. — Eith e r possibilities: it will either roll along a roll corridor formed by 
compressively touching both Spheres 1 and 2 or it will over-roll Sphere 2 and lose contact with 
Sphere 1. To check for compressive contact with both Spheres 1 and 2, the preferred 
implementation checks whether the current sphere would contact Sphere 1 if it were to roll down 
the longitud e longitudinal line of Sphere 2 where the current sphere-has - had j ust contacted 
Sphere 2. If the current sphere would contract contact Sphere 1, then it must be in compressive 
contact with Sphere 1. It is certainly in compressive contact with Sphere 2 because of the 
contact with it. therewith. Therefore, it is in compressive contact with both spheres. If, on the 
other hand, the current sphere would roll away from Sphere 1 while rolling down the longitud e 



longitudinal line of Sphere 2, then it will leave Sphere 1 and is only in compressive contact with 
Sphere 2. 

[00158] Consider the positions of the three spheres: (Xu Yj, Zy) for Sphere 1, (Xi, Y2, Zi) 
for Sphere 2, and (X C9 Y c , Z c ) for the current sphere. The radii of these three spheres are aj, ci2, 
and a c , respectively. The Sphere 2 longitudinal line at which the current sphere is touching 
Sphere 2 is at azimuth 



[00159] If 0 c2 > 7t/2, then the current sphere has rolled into the underside or southern 
hemisphere of Sphere 2. In this case, it will be in compressive contact with the two spheres. If 
9 C 2 < n/2, then we must see if it contacts Sphere 1 if it were to roll down the longitude 



longitudinal line at cp C 2. 

[00160] Let the distance between the current sphere and Sphere 1 be denoted D cI which 
is a function of 0 C 2. Its square is given by 



<p c2 = tm-l(Y c -Y 2 )/X c -X 2 )] 



(16) 



and it touches at latitude 



9 c2 = cos" \{Z C - Z 2 ) la c -a 2 )]. 



(17) 
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D 2 = (X c - X x f + (Y e - y,) 2 + (Z c - Z,) 2 = [(a e + a 2 )sin G c2 cos <p c2 +X 2 - X,f + 
[(a c + a 2 )sm 9 c2 sin <p c2 + Y 2 - 7,] 2 + [(a c + a 2 )cos 0 c2 + Z 2 - Z,] 2 . 



(18) 



If D c] or its square increase as 8 C 2 increases (i.e., as the current sphere rolls down Sphere 2's 
longitud e longitudinal line), then it has lost contact with Sphere 1 . That is, if its derivative with 
respect to 0 C 2 is greater than or equal to zero, it is not in compressive contact with Sphere 1 . If 
the derivative is negative, it is in compressive contact with both spheres. The expression below 
is the derivative to within a positive constant of 2(ac + a2)> so if it is negative, the current sphere 
compressively touches both Spheres 1 and 2: 

[{a c + a 2 )sin 0 c2 cos <p c2 +X 2 - X x ] cos 6 c2 cos (p c2 + [(a c + a 2 )sin 0 c2 sin <p c2 + 
Y 2 - cos 0 c2 sin (p c2 - [(a c + a 2 )cos 9 c2 + Z 2 - Z,] sin 0 c2 . 

[00161] These two-sphere calculations are done in this implementation by Subroutine 

ST2S. 

Compressive Touching Of Three Sphcres Compressive touching of three spheres 
[00162] In the case of compressive touching of three spheres, the current sphere at 
position (X C1> Y C9 Z c ) and with radius a c is in contact with three pack spheres at positions (Xi, Yj, 
Zj) 9 (X 2 , Y 2 , Zi), and (X 3 , Y 3 , Z 3 ) and with radii a u a 2 , and a 3 . Let the current sphere have weight 
W due to gravity or similar force. Assume it is attached to the three spheres at the point of 
contact with them. The preferred implementation solves a force equation to determine whether 
the attachments are compressive or tensile. If an attachment is tensile, the current sphere will 
roll away from it. 

[00163] Let the unit vectors pointing from the current sphere toward Sphere i where i is 
1, 2, or 3, be denoted by-e^r-_e/^ Let the magnitude of the force applied on the current sphere by 
Sphere / be denoted by fi. If fi is positive, the attachment is tensile; if negative, it is compressive. 
The three equations needed to solve for the three yj are given by the vector equation 



//C/ + /a + /a = (0,0,-FT). 
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(20) 



[00164] Because in the implementation we are only interested in whether the fs are 
compressive (ft < 0) or tensile (ft > 0) and are not interested in their magnitudes, the preferred 
implementation simply sets Wto unity and obtains for the / solutions 

// = {e 3x e 2y - e^)/det E 

f 2 = (e lx e 3y - e 3x e Iy )/det E (21) 
fs = {e 2x e ly - e Ix e 2y )/det E 

where det E is the determinant of the 3 x 3 matrix whose columns are the unit vectors-e^- ^, and 
64t e u and e? : 

det E = e lx (e 2y e 3z - e 3y e 2z ) + e 2x (e 3y e h - e Iy e 3z ) + e 3x (e Jy e 2z - e 2y e lz ) . (22) 

[00165] If the determinant is zero, the set of equations has no solution because the three 
placed pack spheres and the current sphere all lie in the same plane. In the unlikely event that 
this should happen in practice, the preferred implementation first sees if another set of three 
spheres gives stability if the current sphere is touching more than three placed pack spheres. If 
there is no stable triplet, then it enacts a subroutine to perturb the current sphere's position 
slightly in a downward direction such that it will not overlap any of the three spheres. If the 
perturbation causes contact to be lost with one or two of the touching spheres, it can proceed 
with FOJB2 or FOJB3. If no contact is lost, then at least the four spheres are no longer coplanar 
and it can redo the processing without a zero determinant. 

[00166] These three-sphere calculations are done in the preferred implementation by 
Subroutine ST3S. 
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Compressive Touching Of One Sphere And The Wal l Compressive touching of 
one sphere and the wall 

[00167] In the case where there is compressive touching of one sphere and the wall, the 
wall cannot be the first object contacted because in this implementation, the current sphere drops 
parallel to it when free falling. Therefore, the only way the current sphere can be in contact with 
both its cylinder wall and one placed pack sphere is if it hits Sphere 1 and subsequently rolls into 
the wall. The current sphere is always in compressive contact when rolling on one sphere. 
Furthermore, because the current sphere must be in contact with Sphere l's northern hemisphere 
when it contacts the wall, it must be in compressive contact with the wall because it just rolled 
into it. Therefore, any time the current sphere is in contact with one sphere and the wall, it is 
always in compressive contact with both and will descend along a roll corridor defined by both. 
In the unlikely event that the current sphere hits the cylinder wall and nothing else at the same 
time it reaches the equator of Sphere 1, it free falls from that point and FOBJ1 is r e call e d. 
recalled. 

Compressive Touching Of Two Spheres And The Wal l Compressive touching: of 
two spheres and the wall 

[00168] The case of compressive touching of two spheres and the wall is a simple 
offshoot of the three-sphere case. The cylinder wall will replace one of the spheres as a contact. 
However, in the absence of friction, the wall can only apply a force to the current sphere directed 
radially inward. 

[00169] Let-^-and-*^ n and n denote the positions of the two pack spheres and iy - and 
rc.denote the position of the current sphere. The relative positions of the two pack spheres to the 
current sphere is-* ^ ~ n r „ - rir = n- r c , i = 1 or 2. These relative vectors point from the current 
sphere center to the z-th pack sphere center. Let-e^-_^_be the unit vector pointing along-^-_r/^ 
The force the current sphere applies to each pack sphere and the wall must sum to its weight W. 
That gives the equation 

F#4^#*r+fa^ - (0,0, m. F ,e,r. + f?e? r + f w e wr = (0. 0. -W). (23) 

41 



[00170] The unit vector^e^H^which points radially outward from the cylinder center, 
is given by 

e^^O^jme}^^ ^ = (xa + rA)Kx} + vh m (24) 

[00171] As before in the three sphere case, the preferred implementation solves for the 
Ji, i = 1, 2, or w, but actually only their sign is needed: positive, negative or zero. Hence Wis set 
to unity. Let be the x component of the vector^— = 1, 2, or w, and similarly for the y 
and z components. Then the solutions for the fx are given by 

// = fay e wx - e 2x e^y)/ det e 

fi = (ej x - e iy O/det e (25) 
fw=(e Ix e 2y - ei y ^/det e 

where 

det e - -e }x e 2z + e iy e 2z e wx + e iz (e 2x - e 2y e wx ). (26) 

If an f h i = 1, 2, or w, is negative, the force between the current sphere and that object is 
compressive and they stay in contact. If it is positive, the force is tensile and the current sphere 
will roll away from that object. In the unlikely event that it is zero, the implementation considers 
it tensile and has the current sphere roll away from the object, thereby trying to find a point of 
stability at lower altitude. If all three f are negative, the sphere is in stable contact and is placed 
at its current position. 

[00172] These stability calculations for two spheres and a wall are performed in the 
preferred implementation by Subroutine ST2S1W. 

Detail of Subroutine FOBJI Detail of Subroutine FOBJ1 

[00173] Returning to the subroutine FOBJ1, it finds the first object the descending 
current sphere encounters as it is "dropped" from the position selected by subroutine 
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DROPSITE. There are three possible objects the current or selected sphere can encounter: (1) an 
already-placed pack sphere, (2) a water level, or (3) the catch net. A fourth object that may be 
encountered later is the cylinder wall, but because the current sphere is assumed to drop parallel 
to the cylinder wall, it cannot contact the wall until it rolls off some placed pack sphere. 

The Catch Nc t The catch net 

[00174] In accordance with the method aspect of the invention, the selected particle 
placement includes defining a catch net a catch net representative of buoyancy o f a portion of 
the placed particles a portion of the placed particles and positioning the catch net within the 
space based upon the placement of the portion of the placed particles. The selected particle is 
placed in non-overlapping relation with respect to the catch net. Preferred versions of the 
method include defining a pack surface for the placed particle, and positioning the catch net 
relative to the pack surface. The pack surface may comprise a top layer of the placed particles 
within the space, or preferably, within each of the subspaces. The particle surface optionally 
may correspond to an average of the south poles or south pole positions of the top layer of placed 
particles. The catch net may be positioned at a distance from the pack surface based upon a 
selected particle radius, such as an average particle radius of the top layer or another subset of 
the placed particles at the pack surface. 

[00175] In the way of background, when simulating a particle pack, particularly one 
with a wide particle size distribution, if the volume fraction of large particles is much larger than 
that of small particles, the small particles can rattle down through the pack by falling through the 
interstitial spaces of the large particles. In many particulate materials, there are buoyancy or 
suspension considerations that would keep the small particles more or less well distributed 
throughout the pack. Examples of such buoyancy or suspension factors may include the nature 
or properties of the binder or other inter-particle medium or media, the interaction of the particles 
with the inter-particle medium or media, interactions among the particles themselves (such as 
electrical or magnetic interactions), interactions with e xt e rnals external forces or fields that tend 
to disperse or suspend the particles, and similar phenomena. 

[00176] In accordance with this aspect of the present invention, this buoyancy or 
suspension phenomena is simulated using a "catch net." The catch net comprises an assumed 
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surface or construct to catch the current or selected particle as it is placed and to prevent that 
particle from passing a defined location or level which theoretically represents space occupied by 
smaller ones of the previously placed particles. 

[00177] In the presently preferred implementation, a catch net is incorporated into the 
space, below which the south pole of the falling particles cannot descend. When a particle's 
south pole hits the catch net belonging to its cylinder, it is placed in the pack at that point, even if 
it is not touching any other particle. Preferably there is a separate catch net for each particle 
submode, and thus for each subspace or cylinder. Each catch net preferably covers the entire 
cross section of that submode's cylinder. If the current sphere's south pole hits that catch net, it 
is placed there. It does not see or respond to any other submode's catch net. The altitude of the 
catch net preferably is determined dynamically, that is, it can depend on which placed particles 
are being contacted. 

[00178] The use of a catch net is not necessarily limited to a specific type or size of 
particle, or a particular size or geometry of the space or subspaces. In accordance with the 
presently preferred implementation, however, as noted, the particles comprise spheres spheres 
and the characteristic dimension of each of the particles comprises a radiu s , a radius. The 
central string comprises a lin e , pr e f e rably the z axi s , a line, preferably the z axis, within the 
space. The space comprises a cylinder, a cylinder, and the subspaces comprise concentric 
cylinders concentric cylinders disposed about the central string, each of the cylinders defining 
one of the subspaces and having a radius a radius that corresponds to the radius of the particles 
in the category corresponding to that cylinder. 

[00179] In accordance with the presently preferred implementation, the method further 
includes defining a pack surface a pack surface for the placed particles, and the catch net 
positioning comprises positioning the catch net relative to the pack surface. Also in accordance 
with this implementation, the catch net positioning comprises positioning the catch net— a 
di s tance from the pack s urface a distance from the pack surface based upon a s elected 
particle radiu s , a selected particle radius. 

[00180] The placed particles may be assumed to comprise a top layer, a top layer, and 
the pack surface optionally but preferably comprises an average o f the particle location s the 
particle locations of the placed particles in the top layer. In this implementation, the particle 
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surface corresponds to the south poles of the top layer placed particles and more specifically 
corresponds to an average of the south pole po s ition s an average of the south pole positions of 
the top layer placed particles. 

[00181] In the presently preferred implementation, the catch net extends across the cross 
sectional area of the space, across the cross sectional areas of each of the subspaces. 

[00182] In this preferred implementation, the catch net comprises N s ubn e t s , N subnets, 
one of the subnets corresponding to each of the subspaces. Each of the concentric cylinders 1,2, 
3, .. N I 2, 3, ... N has a circular cross section perpendicular to the central string. In the 
preferred implementation, the catch net comprises a plurality N of subnets, or individual catch 
nets, one for each of the cylinders. Each of the subnets preferably extends over the cross 
sectional area of the corresponding subspace or cylinder. The position or level of the catch net or 
subnet for each of these subspaces can and typically will differ from one another, although in 
some instances they may be at the same or nearly the same level. 

[00183] Preferably, the positioning of the catch net positioning comprises setting the 
catch net or each of the subnets at a s e lected a selected distance from the top surface. This 
selected distance may be a user-defined quantity, or it may be dictated by the material being 
simulated or other factors. The catch net preferably is positioned at an altitude along the z-axis 
that is a function of the particle radius of particles corresponding to the mode or submode of the 
cylinder in which the catch net or subnet resides. The catch net or subnet pr e f e rably but 
optionally preferably, but optionally, is some multiple of the particle radius, preferably one or 
two ratio, below the top layer or pack surface. Adjustments or offsets, however, may be used. 

[00184] When the pack is first being formed so that not enough particles have yet been 
placed to adequately define a pack surface, the catch net may be defined as the bottom of the 
cylinder into which the particles are being dropped. In some instances, however, particularly 
where the specific application involves simulating a very large three-dimensional pack, it is often 
advantageous or desirable to avoid the edge effects of the flat bottom of the cylinder. 
Accordingly, in some instances the positioning of the catch net ma y comprises comprise spacing 
the catch net from the base surface of the space. The spacing from the base surface, for example, 
may simulate the positioning of the catch net fo r a top layer of the placed particle s , a top layer 
of the placed particles, for example, the positioning of the catch net off of the base of the space 
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may be done in a way that the bottom layer of placed particles simulate the shape, configuration 
and appearance of the top layer of particles after the particle pack has been constructed to include 
multiple layers. 

[00185] This may be implemented by assigning as initial catch net po s ition for a kth 
one of the s ub s poce s an initial catch net position for a k-th one of the subspaces Z\ n \ t (k), 
assigning as th e charact e ri s tic dim e n s ion of the particle s of a kth on e of the particle 
categories the characteristic dimension of the particles of a A>th one of the particle categories a^ 
assigning as the characteristic dimension o f a s mall one of th e particle s a ^ of th e particle s , _ a 
small one of the particles amin of the particles and assigning as the characteristic dimension of-a 
large one of the particle s a w ^- a large one of the particles a,^. Below a threshold level _a 
threshold level the catch net positioning further comprises positioning the catch net fo r a kth one 
of the particle categories at a catch n e t position a £-th one of the particle categories at a catch 
net position Zn Pf (k) within a kth on e of the s ubspaces a £-th one of the subspaces determined by 

Z net (k) = Z init + Hra k a min /a max (27} 

where r represents a weighting co e ffici e nt a weighting coefficient and H represents a s witching 
coefficient, a switching coefficient. The weighting coefficient pref e rably but optionally 
preferably, but optionally, is assigned a random number. Below the threshold value the 
switching coefficient preferably is assigned a value of one, and above the threshold value the 
switching coefficient preferably is assigned a value of zero. 

[00186] In the presently preferred implementation, the switching coefficient mechanism 
is implemented by a switch called ROUGHNET which, when true, invokes the random number 
generator to slightly change the altitude of the catch net so that the first set or layer of placed 
particles do not all have the same altitude. If ROUGHNET is false, the south poles of the first 
set of placed particles do have the same altitude because they hit the same catch net. 

[00187] Let Z init be the initial altitude above which the pack will be built. It is specified 
by the user in the preferred implementation but its default preferably is zero. The initial catch 
net altitude for the k-th mode is 
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Znet(k) = Zinit + H r a k a min /a max 



mm 



where r is a random number in [0,1], \0. 11 a m i» is the smallest particle radius of all submodes 
belonging to the pack, a max is the largest, and H is unity if the ROUGHNET switch is true and 
zero if it is false. This expression for the catch net altitude holds until enough particles have 
been placed to establish a pack surface, i.e., until m (which is greater than or equal to m min ) is 
reached. No part of any particle is permitted to descend below Z init in this implementation. 

[00188] The catch net preferably is positioned relative to the pack surface. The pack 
surface can be determined in a number of way s . It may be determined, for example, ba s ed 
on th e large particle s in the top layer, or another s ub s et of the placed particles, e.g., such as 
the la s t particles placed, can be determined in a number of ways. It may be determined, for 
example, based on the large particles in the top layer, or another subset of the placed particles, 
e.g., such as the last particles placed. 

[00189] The particle surface may be ascertained in the following manner. For each 
particle category k particle category k of the placed particles in the top layer, the program 
defines a particle radius o r a particle radius ai for the plac e d particle s i placed particles i of 
that category k. For the cylinder k corresponding to the particle category k, the programs 
program assigns a cylinder radius W k- a cylinder radius Ww. The program also assigns a top 
layer particle number m(k) and determines values for m(k) by evaluating 



m(k)-\ m(k) 

I a]<WZ< Zaj,k = \,2 N .... 

/=1 /=! {Zy ) 

Submodel)** Submode(i)Zk 



where N is the number of particle categories, e.g., the number of submodes. The program then 
determines the pack surface location using 



1 m 

S^ — TiZi-ad (30) 
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where S represents the pack surface and Zj represents the position of a center o f a cent e r one of 
the placed s ph e r es , a center one of the placed spheres. 

[00190] In accordance with the preferred implementation, the surface altitude S of the 
pack is defined as the average south pole altitude of the last m placed particles where m is a 
positive integer that will be discussed below. This surface altitude is given by 



1 m 
m '=i 



where Z, is the altitude of the center of Sphere i which is one of the last m placed particles and at 
is its radius. Instead of using a sinRle4th - z-th p article radius in Equival e nc e . 27, Eq. 31, one 
may generalize the expression by permitting (Zi - aO to be written (Zj - f(aO) where f(aO is a 
function of particle radius a;. The preferred value for f(aO = a i? but f(aO may assume other forms, 
for example, such as f(ai) = na; where n is a positive real or integer number. 

[00191] Although the "surface" of a random particle pack may-be-a- be defined in a 
number of ways, the surface preferably is defined only after the integer m is large enough that the 
falling spheres have covered the cylinder area. The k-th cylinder with radius W k has area n W k 2 . 
Let m be defined as the number of most recently placed spheres whose submodes have cylinder 
radii less than or equal to W k and whose summed cross-sectional area just exceeds or equals the 
&-th cylinder area. For the A:-th submode, m is taken from the following condition. Let m(k) be 
the number m for the &-th submode. It is found from 



m{k)-\ m(k) 

Z^<»l> S g ai >. 4 = 1,2 N 

Submode(i)Zk Sub mod e(i)<.k 



where TV is the total number of submodes. The submode k for which this condition is satisfied 
first, that is, has the lowest m(k) becomes the m for the pack: 

m - MTN[m(k)l (29)03) 
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[00192] In the preferred implementation there is a minimum number m min for those 
unused packs that defeat the reasoning leading to the value for m given in Eqs (28) and (29), In 
th e cas e s th e Eqs. 32 and 33. The default value of m m j n to 12 is may be assigned to a value of 
12 in the preferred implementation. 

[00193] Once a surface is established, its altitude S is determined by using Eq. (27). 
Eq. 31. As the current sphere descends, if it hits a pack sphere much larger than itself, and if the 
TUNNEL switch is true, it will tunnel through the large sphere to check for possible cavities 
under the large particle. Therefore the catch net should be placed quite low so the current sphere 
can investigate possible cavities. If TUNNEL is false, the catch net does not need to be that low. 

[00194] If the current sphere hits a sphere somewhat bigger than itself, then it is unlikely 
to tunnel because the surface sphere is not large enough to accommodate it in a cavity. Then the 
catch net need not be as deep. 

[00195] Finally, if the current sphere hits a sphere smaller than itself, then it cannot 
penetrate too deeply into the pack so the catch net can be rather shallow. 

[00196] Letting the subscript i denote the first pack sphere contacted, the catch net 
criteria are quantitatively summarized as follows: 

[00197] If aj/a c < 1, then the catch net is shallow: 

Z net {k) = S-a h (£±a)(34a) 
[00198] If 1 < aja c < a Xi where a x preferably is V6 + 2, then the catch net is moderate: 

Zneik) = S-2a c . Q4b ¥34b) 
[00199] If ai/cic > a x , where a x preferably is V6 + 2, then the catch net is deep: 

Z net {k) = -2a c - at . B±e ¥34c) 

[00200] The value of a x represents the sizes required for a current sphere to fit into a 
cavity formed by larger spheres, which preferably assumes a value of (V6 + 2). If the centers of 
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four identical larger spheres form a regular tetrahedron, and a smaller sphere just fits into the 
interstitial space of the tetrahedron, then the ratio of the smaller sphere to larger sphere diameter 
is V6 + 2 in preferred implementation. 

[00201] Further in accordance with the preferred implementation, the code monitors the 
selected or current particle placement to identify the occurrence of the following sequence of 
events: 

[00202] 1. The catch net is set at an altitude when the current sphere hits the 

z-th pack sphere. 

[00203] 2. Then the current sphere rolls off the i-th sphere and loses contact 

with it. 

[00204] 3. When it hits another pack sphere, the catch net is set above the 
current altitude of the current sphere's south pole. 

[00205] If this series of events is detected, the current sphere is raised (its position is 
adjusted upward) until its south pole is at the new catch net altitude, and it is then placed there. 

The Water Levc l The water level 

[00206] In accordance with another aspect of the invention, the placement of the 
particles comprises defining a water level a water level representative o f a level of a portion of 
th e placed particle s that are s maller than the s el e ct e d particle a level of a portion of the 
placed particles that are smaller than the selected particle and represent a s urface a surface of 
the smaller placed particles, and positioning the water level within the space based upon the 
smaller particle surface. The current or selected particle4heft is then p laced in non-overlapping 
relation with respect to the water level. 

[00207] The water level serves a purpos e different purpose than the catch net. Recall 
that the particle pack according to the preferred implementation is constructed using a series of 
concentric cylinders, one cylinder defined for each particle category or submode. Recall further 
that the centers of particles in this implementation cannot pass outside their own corresponding 
and constraining cylinder into which they are dropped. Large cylinders are used for large 
particles and small cylinders for small particles, unless the user for some reason designates 
otherwise. The cylinder construct facilitates or permits use of the reduced dimension approach 
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to particle pack construction. Its purpose is to prevent the program code from having to fill large 
cylinder volumes with potentially very large numbers of small particles. 

[00208] This approach, however, creates the problem that the large cylinders will 
contain empty space into which large spheres can roll during simulated particle placement where, 
in a corresponding physical particle pack such as the one being simulated, this would not occur. 
In the actual pack, the smaller particles would be included even in the regions comprising the 
larger cylinders, and these smaller particles would be encountered to prevent such rolling. 

[00209] The water level simulates this presence in the larger cylinder areas of smaller 
particles that, according to the presently preferred implementation, have not been placed in the 
larger cylinders. 

[00210] There are a number of methods or approaches for selecting or setting the water 
level. According to one such approach, the water level positioning comprises determining-aa 
average location an average location of the particle locations along the central string of the 
particles o f th e portion of placed particles the portion of placed particles that are smaller than 
the current or selected particle, and positioning the water level at the average location. The water 
level may be set at the average altitude o f the small particle " s urface" the small particle 
"surface" as if the outside cylinders were filled with small particles. Then when a large particle 
or sphere descends onto the water level, the particle, and preferably its south pole, is stopped 
there as it would be if that volume had been filled with the smaller particles. 

[00211] In the presently preferred implementation, wherein the plurality of subspaces 
(here concentric cylindersVare- is used, the water level comprises a plurality of subspace water 
levels, preferably one subspace water level for each subspace or cylinder. These water level 
positions may and typically will differ from one another as the particle pack is constructed. The 
water level positioning optionally but pr e f e rably optionally, but preferably, comprises assigning 
a s ubspace water l e v e l po s ition to a portion a subspace water level position to a portion of the 
subspace water levels. The portion of the subspace water levels preferably correspond to one 
water level for each cylinder except the first cylinder. In the preferred implementation, all of the 
cylinders except the first or central one, containing the smallest particles, has a subspace water 
level and corresponding subspace water level position. 
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[00212] In the presently preferred implementations, the water level positioning 
comprises using an offset an offset to position the water level, o r an off se t an offset relative to 
the average location average location of the particle locations. Preferably, the water level 
positioning comprises using-on - an offset for each of the s ubspoce water level position s . 
subspace water level positions. The water level positioning similarly may comprise using-a» 
offset - an offset for each of the s ubspac e water level po s ition s , subspace water level positions. 

[00213] In the presently preferred implementation, the water level altitude for the i-th 
cylinder (i> 1) is defined as the average exposed north pole altitude of all particles with radius 
less than ai. The user may raise or lower this level by adding a user-defined offset, but the 
default for this offset in this implementation is zero. This water level exists between the /-th 
cylinder wall of radius Wi and the next smaller cylinder wall of radius W^. Thus the spaces 
between different cylinder walls have different water levels. The cylinder corresponding to the 
smallest particle's submode has no water level, only a catch net. 

[00214] If a larger sphere is descending onto the pack near the center cylinder, it may hit 
the center cylinder's pack of small particles and tend to roil off its edge and fall to the floor of 
the pack or fall until it hits another large sphere. Thus the pack would segregate with larger 
particles outside smaller particles. With the water level in place, however, if the larger sphere 
rolls off the pack in the central cylinder, it simply "floats" on the water level which is at the same 
level as the "surface" of the central cylinder pack. Thus, any time a descending sphere hits a 
water level, it is placed where it is, just as was done with the catch net. Ordinarily, water levels 
will be higher than catch nets, although this is not always true. 

Search for the First Object Encountcre d Search for the first object encountered 
[00215] A free-falling current sphere can encounter the following situations. It might hit 
a single pack sphere. If so, FOJB2 is called. It can hit two pack spheres simultaneously. If so, 
FOBJ3 is called. It can hit three or more pack spheres simultaneously. If so STABLE is called 
to see if it is in stable contact with any triplet of these spheres. Hitting two or more pack spheres 
simultaneously is highly improbable in a double precision random code but not so improbable if 
the pack is semi-ordered, and hence these scenarios are included in the preferred implementation. 
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[00216] In addition, the south pole of the current sphere may hit the catch net or a water 
level before it hits any other sphere. In either case, it is placed where it is. 

[00217] A pack sphere, e.g., the y'-th sphere, is a candidate for being the first object 
encountered by the descending current sphere if their two radii overlap, i.e. 9 if 

(x c - xj) 2 + (y c - Yj) 2 < (oc + a/ mm 

If there were no other obstacles in the way, the current sphere would strike Sphere j at an altitude 
of 

Z c (j) = Zj + [(a c + ajf • (X c -Xj) 2 - (Y c - Yj) 2 f 2 . mX36) 

[00218] Once the pack has grown to many spheres, efficiency can be greatly enhanced 
to the extent that as many spheres as possible be excluded from the detailed considerations of the 
above two equations, while not unduly adversely affecting the accuracy of the simulation in 
representing the physical particle pack being simulated. For example, if | x c - xj | or | y c - yj \ 
alone were greater than (a c + aj) 9 then there is no need to check the square of their sums as 
above. 

[00219] The processing demands are reduced in the preferred implementation because 
all spheres whose north poles have submerged beneath the catch net are off limits and need not 
be included. Simply flagging them by putting a minus sign in front of their mode or submode 
number indicates that no calculations are to be done to see if they are in range of the descending 
current sphere. Furthermore, in this implementation the program keeps track of the lowest index 
number for any sphere whose north pole protrudes above that mode's catch net so that the 
eventually large list of particles beneath the catch net are not searched at all. 

[00220] Submergence of a placed pack sphere beneath the water level, however, is not 
necessarily grounds for dismissing it from consideration. The current sphere may roll into a 
cylinder for which the water level is different and the pack sphere may again become a candidate 
for contact. 
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[00221] If the current sphere is dropped far enough outside a set of cylinders that there is 
no chance for a first contact with the spheres belonging to those cylinders, then in the preferred 
implementation-thev - such sphere(s) are dismissed from further consideration. If the submode 
number of a sphere belongs to an out-of-range cylinder, no further processing is required in the 
preferred implementation. Let W k be the cylinder wall radius for the &-th particle submode. If 

{W k + a k ) 2 <X c 2 + Y C 2 (34)(37) 

then no k-th mode sphere is a candidate for first contact with the descending current sphere in 
this implementation. 

[00222] In summary then, FOBJ1 searches the list of all placed pack spheres above the 
catch net for the submode to which the current sphere belongs, excluding all placed spheres that 
may be in interior cylinders that the current sphere cannot contact or overlap. The possibility of 
overlap in the x directions is checked, followed by the possibility of overlap in the y directions, 
followed by a direct calculation of overlap. The altitude z c (j) at which the current sphere would 
contact Sphere j is stored and that sphere is rejected if other placed pack spheres have a higher 
contact altitude. Again, if no pack spheres make contact before the current sphere reaches the 
catch net or water level, it stops and is placed there. 

Detail of Subroutine F0BJ 2 Detdi\ of Subroutine FOBJ2 

[00223] Subroutine FOBJ2 finds the second object to be contacted by the descending 
current sphere. FOBJ2 would only have been called if the current sphere had first contacted a 
single pack sphere that has been denoted Sphere 1. If it had contacted the catch net or water 
level, the current sphere would have been placed there and FOBJ2 would never have been called. 
If it had contacted more than one sphere simultaneously, FOBJ3 and/or STABLE would have 
been called. FOBJ2 is called only if a single pack sphere's northern hemisphere has been 
touched. 

[00224] When the current sphere touches the northern hemisphere of Sphere 1, it begins 
to roll down the path of steepest descent which is Sphere P s longitudinal line pa ss ing through 
the point of contact. Sphere l's longitudinal line passing through the point of contact. The 
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volume swept out as the current sphere descends along this path toward the equator is called the 
roll corridor. The end of the roll corridor is either where the current sphere's center reaches the 
equatorial plane of the Sphere 1, reaches its cylinder wall, or reaches the catch net or water level, 
whichever is higher. The preferred implementation identifies all other pack spheres that protrude 
into the roll corridor. 

[00225] Processing can be somewhat simplified if the cylinder wall is assumed to be the 
boundary that the center of the current sphere cannot pass, rather than the boundary through 
which its leading edge equator cannot pass. 

[00226] If the center of the current sphere intersects the cylinder wall during descent, 
then that point is the end of the roll corridor because it will at least have found the wall as its 
second object if no other pack sphere has been contacted first. 

[00227] The beginning of the roll corridor is then the place of initial contact between the 
current sphere and Sphere 1. The current sphere's polar angle at that initial contact is given by 



[00228] The end of the roll corridor 9 c f will either be at the equatorial plane of Sphere 1 
or where the center of the current sphere intersects the catch net, the water level, or the cylinder 
wall. If it is at the equatorial plane, then 



G ci = cos" [(Z c - Z } )/(a c + ai)\ 



0 c f = 7t/2. 



W£39) 



If the altitude of the catch net Z net >Zj - a c , then 



6 C f = cos' [(Z net + a c - Zj)/(a c + a/)]. 



mm 



If the altitude of the water level Z w >Zj - a C9 then 



0cf = cos' } [(Z w + a c - Zj)/(a c + a/)]. 



mm 
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[00229] The last possibility is that the roll corridor ends with an intersection of the 
cylinder wall. Let cp c be the azimuthal position at which the current sphere has touched Sphere 1 . 
It is given by 



<p c = tan to - Yj)l{X c - X } )1 (49)(42) 



[00230] The radial position of the current sphere as it descends Sphere 1 as a function of 
0 C is given by 

{[Xj + (a c + ai) sin6 c coscp c ] 2 + [7/ + (a c + ai) sin6 c sin0 c ] 2 } 1/2 . (4©)(43) 

[00231] If this radial position is equal to the cylinder radius W for some value of 0 C 
between 0 and n/2, then the roll corridor intersects the cylinder wall. Setting this expression 
equal to W, squaring both sides, and solving the quadratic equation for sin0 c gives 

. _ -(X x cos (p c + y, sin <p c ) ± ^{X x cos <p c + Y x sin <p c ) 2 + W 2 - X x 2 - Y? 



[00232] Because 0 c f is a polar angle, it lies in the range [0,7t]. [0, tc]. Therefore, sin0 cf 
can never be negative. If both right hand right-hand side solutions o f Eq. ( 4 1) Eq. 44 are 
greater than unity, there is no solution and the current sphere does not intersect the cylinder wall. 
If both are less than or equal to unity, then we choose the larger one. 

[00233] The presently preferred implementation conducts appropriate processing to 
evaluate each of these possible circumstances. Before evaluating Eq. ( 4 1), Eq. 44, the program 
code according to this implementation adds the global radial position of Sphere 1 to the sum of 
the radii of Sphere 1 and the current sphere. If this sum is less than or equal to W, then the 
sphere cannot hit the cylinder wall in any case so the more involved evaluation above is 
unnecessary and preferably is not performed. 
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[00234] The end of the roll corridor is the minimum 0 c f o f Eqs. (36, 37, 38, and 4 1). 
Eqs. 39, 40, 41, and 44. If the current sphere encounters the catch net or water level, it is placed 
there. 

[00235] The preferred implementation also checks for all pack spheres whose volumes 
intrude upon the roll corridor swept out by the current sphere which covers the arc from 0 C j to 0 c f 
and azimuth (p c on Sphere 1. If there is a 0 C with smaller value than the end of the roll corridor 
9 C f at which the current sphere encounters one or more pack spheres, then the sphere or spheres 
touched are noted. Program control is then passed to subroutines that deal with simultaneous 
contact of two or more spheres, as described herein. 

[00236] If no pack spheres overlap the roll corridor, then the current sphere is placed if it 
hits the catch net or water level. If it hits neither of these, then it rolls to the equator and drops 
from there. At that point it is a free falling sphere again so FOBJ1 is r e call e d, recalled. 

[00237] We now turn our attention to finding which placed pack spheres overlap the roll 
corridor and at what point the current sphere would strike them as it descends from 9 C i. The term 
"global coordinate system" is used herein when discussing a single coordinate system in which 
the entire pack is defined. Capital letters are used when referring to a global coordinate system, 

e.g., Ri =-(Xx¥xZi}- (Xu Y,\ Zj) for the Cartesian position of the z-th sphere. 

[00238] A local coordinate system, on the other hand, is one whose origin is at the center 
of one of the spheres and whose axes are parallel to the global axes. We use small letters to 
define it, e.g., r; ^H^r^rh fc, y/ , zj) for a position relative to the center of the z-th sphere. 

[00239] Consider the y'-th pack sphere. Its center is at global position-J?/- J^_= (Xj, Yj, Zj) 
and its radius is oj. In the global coordinate system, the position of the current sphere as it 
descends along Sphere / is 

R c = + (a c + a s ) (sin 0 C cos <p c i + sin 6 C sin (p j + cos 0Jc) • (42){45) 

[00240] The current sphere rolls into Sphere j if there is a solution 9 C that lies in the 
interval 9 C i <9 C < 9 C f to the equation 
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R c - Rj 



mm 



[00241] If this equation can be satisfied at all, it will have two solutions: one where the 
current sphere touches Sphere j from above and another where the current sphere touches Sphere 
j from below. (If it just nicks Sphere j, these two solutions can be at the same 8 C , but if Sphere j 
is just nicked, the preferred implementation considers that a miss.) We are only interested in the 
top or smaller 9 C solution, and then only if the top 6 C lies in the roll corridor range. 

[00242] Explicitly, the equation to be solved is 

[X t -Xj + (a c + ai) sin0 c cos<p c ] 2 + [7/ - Yj + (a c + ai) sin0 c sincp c ] 2 + 

[Zt - Zj + (a c + ad cose c ] 2 = (a c + aj) 2 (44)(47) 



which, defining Xy = X x • • Xj, Yy = 7, - Yj, and Zy = Z,- - Zj, can be rewritten 



(Xy cos(p c + Yy sincpc) sin0 c + Zy cos9 c = [(aj - a t ){aj + a { + 2a c ) - Ry 2 ]/[2(a c + ai)] (4§){48) 



where 



Ri=Xj + Yy 2 + Zi, 2 . (46)(49) 



This equation has the form 



^ sin0 c + B cos0 c = C (47)(50) 



where 



^4 = cos(p c + 7^ sincpc 
5 = Zy 

C = [(aj - ad(aj + a, + 2a c ) - Ry]l[2(a c + a,)L 



(49)(52) 

mm 
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Dividing through by ^l(A 2 + B 2 ), we obtain 



cosa sin6 c + sina cos9 c = CH(A 2 + B 2 ) f§4¥54) 



where 



cosa = AN(A 2 + B 2 ) f52¥55) 



and 



sina = BN(A 2 + B 2 ). m*56) 
[00243] A common trigonometric identity turns Eq. (51) Eq. 54 into 

sin (0 C + a) = CN(A 2 + B 2 ) (54){57) 



which has solution 



0 C = sin 1 [CN(A 2 + B 2 )] - a 



where 



a = tan-'[^/5]. W£59J 

If | CN(A 2 + B 2 ) | > 1, then there is no solution and the y'-th sphere does not overlap the roll 
corridor. 

[00244] The FORTRAN function ATAN2llff) ATAN2 (A, B) or its equivalent in other 
programming languages uses information on the signs of A and B to determine in which quadrant 
a lies so it is uniquely defined. The arcsine, however, has two solutions. One solution is the 
principal value which lies in \ tc/2,7t/21 r-rc/2, n/2] and is found by the standard FORTRAN 
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function ASIN or its equivalent. The other solution is n minus this principal value. These two 
solutions comprise the two solutions for the current sphere contacting Sphere j on either side 
along the longitud e longitudinal line defined by azimuthal angle cp c . Denoting the principal 
value by capitalizing the arcsine, the two solutions are 



e c i = Sin _1 [C/V04 2 + ^)]-a 



(S7a)(60a) 



and 



9 C 2 = 7i - Sin" 



[CN(A 2 + B 2 )] 



a. 



[00245] If either of these 0 lie outside \02n\ [U 27iL then 2n must be added or 
subtracted enough times to bring the solutions back into r0,27i). TO, 2n]. If neither of these 
solutions lie in the open interval (9 C i, 0 C f), then the current sphere does not contact Sphere j. If 
one solution does lie in (0 C i, 0 c f), that is the position at which the second contact is made. If both 
do, the lesser of the two solutions is the position of second contact. 

[00246] Finally, whichever pack sphere has the smallest 0 C is Sphere 2. That point of 
contact is denoted 0 c j. If two or more pack spheres are contacted simultaneously, then all their 
index numbers are returned by FOBJ2. If the cylinder wall and one or more pack spheres are 
contacted simultaneously, their index numbers are returned and the flag indicating a wall contact 
is activated. 

[00247] As the preferred implementation searches through the pack spheres to find 
which is touched first, there are several efficiency criteria that can be used to eliminate most 
spheres from the detailed consideration of the above equations. 

[00248] First, if the roll corridor is far enough away from the global Z axis, then the 
current sphere may not be able to approach some inner cylinders closely enough to touch their 
spheres. Therefore, those pack spheres whose mode number indicates they belong to the inner 
cylinders need not be considered further. The preferred implementation can determine which 
cylinder's spheres can never intersect the current sphere's roll corridor? corridor. 
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[00249] Because the current sphere rolls down a longitud e longitudinal line on Sphere i, 
when viewed from above, its projection is a straight line segment on the XY plane with beginning 
point (X ch Y cl ) and ending point (X c f 9 Y c j) where 

X ci = Xi + (a c + at) sin9 ci cos(p c f§ga )(61a) 

Y C i = Yi + (a c + ai) sin6 ci coscp c (58b){61b) 

X d f=Xi + (a c + a/) sin0 cf cos(p c (§«e)(61c) 

7 C / = 7/ + (a c + a,) sin0 cf coscp c . (#^ (61d) 

[00250] To find whether this line segment comes close enough to the l-th cylinder to 
touch any of its spheres, the preferred implementation turns the current sphere's position into a 
pair of parametric equations: 

X c - X ci + (X cf -X ci )s (S9a)(62a) 
Y c = Y ci + (Y cf - Y ci )s f59 W62b) 

where s is a dimensionless parameter that is zero at the beginning of the roll corridor and unity at 
the end. In order for the current sphere to have any chance of touching an /-th mode sphere, it 
must come closer than the sum of its radius plus an l-th sphere's radius from the /-th cylinder 
wall which has radius W\. In other words, some portion of the roll corridor must come within 
W\ + at +a c of the global Z axis. This condition is expressed by 



X c 2 + Y c 2 = (Wj +ai + a c f m)(63) 

for a value s somewhere in the open interval (0,1). If Eq. (60) TO, 11 If Eq. 63 is never satisfied, 
then the current sphere cannot touch any of the /-th mode spheres. Explicitly written, Eq. (60) 
Eq. 63 is a quadratic equation in s: 

As 2 + 2Bs + C=0 f6t¥64) 
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where 



a = (x cf -x ci ) 2 + (Y cf - y ci ) 2 mm 

B = X ci (X cf -X ci ) + Y ci {Y cf - Y cl ) mm 

C = Xj + 7 C/ 2 - (Wi + a/ + a c ) 2 . f64¥67) 

[00251] The solutions are 



~B± JB 2 -AC 

s = \ (6SX6S) 

A 



If B 2 -AC is negative, there is no solution. If it is zero, there is only one real solution which 
means the roll corridor just grazes the l-th cylinder, in which case it is considered a miss. If it is 
positive, there are two real solutions. The preferred implementation solves for both of them, and 
checks to determine whether either one lies in the open interval-ffi4V [0, 1]. If one or both lie in 
this interval, the /-th mode spheres are candidates for touching the descending current sphere. 

[00252] To find which cylinders are inaccessible to the current sphere, the preferred 
implementation starts with the innermost cylinder (7 = 1) and continues outwardly until it 
reaches a cylinder that is big enough that it is accessible. All larger cylinders also will be 
accessible. 

[00253] There also is a simple altitude criteria to delimit detailed consideration of pack 
spheres as candidates for touching the current sphere. The global altitudes of the beginning and 
end of the roll corridor are 

Z ci = Zi + (a c + ai) cos9 ci <£6a ¥69a) 
Z cf = Zi + (a c + a t ) cos9 c f , f66b ¥69b) 

respectively. Therefore, if a y'-th pack sphere lies below the end of the roll corridor by the 
amount 
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Zj<Z c f- dc'dj 



(67a)(70a) 



or lies above the top of the roll corridor by the amount 

Zj > Z ci + a c + aj , (63b)(70b) 

then it is not a candidate for touching the current sphere. 

[00254] Likewise, spheres lying beyond the roll corridor in the x and y directions need 
not be further considered, i.e., if 

Xj<X cf -a c -aj (67e)(70c) 

Yj<Ycf-a c -aj (674){70d) 

or 

Xj>Xti + a c + aj (67e)(70e) 

Yj>Y ci + a c + aj, <67*){7O0 

then the y-th sphere lies beyond the roll corridor in the horizontal plane and is not a candidate for 
touching the current sphere. 

Detail of Subroutine FOA/3 Detail of Subroutine FOBJ3 

[00255] Subroutine FOBJ3 finds the third object to be contacted by the descending 
current sphere when it is already in contact with two and only two other objects. There are only 
two possibilities for the two initial contacts: (1) both are placed pack spheres or (2) one is a 
placed pack sphere and the other is the cylinder wall for the submode of the current sphere. Each 
of these possibilities will now be addressed in turn. 
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Touching two sphercs Touchins two spheres 

[00256] It is presumed at the outset of this discussion that the two touching spheres, i 
and j\ are in stable (compressive) contact, subroutine STABLE already having been called to 
remove them if they are not. A roll corridor follows a circle defined where the current sphere is 
in contact with both Spheres / and j. Of course, the current sphere cannot touch both spheres 
unless the condition r 7/ < a, + aj + 2a c is satisfied where ry, is the center-to-center distance 
between Spheres / and j. 

[00257] The current sphere's descent along this roll corridor will be defined by a single 
variable q> c '. This variable is the azimuthal angle of a local tilted coordinate system whose origin 
is in the center of Sphere i 9 with z axis tilted so it passes through the center of Sphere j, and with 
x axis swung upward so that it lies in a vertical plane. In that coordinate system (which is 
denoted herein by primes), q> c ' alone describes the current sphere's position along the roll 
corridor. In the unlikely event that the current sphere finds itself perfectly balanced on Spheres i 
and j so that cp c ' = 0, then a random number generator determines in which direction the current 
sphere will roll - toward negative cp c ' or toward positive (p c '. 

[00258] To find the coordinates of a pack sphere, say the &-th sphere, in this primed 
coordinate system, the preferred implementation first translates the global sphere positions to the 
unprimed local coordinate system whose origin is at the center of Sphere / but whose axes are 
parallel to the global axes: 

r k i=R k -Ri- mem 

[00259] Capital letters are used herein to denote coordinates in the global coordinate 
system, small letters in the unprimed local system (local to Sphere /), and primed small letters in 
the local tilted coordinate system. 

[00260] The Euler matrix used to rotate the unprimed into the primed coordinate system 
is given by 
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- cos 6ji cos <pji - cos 6ji sin (p^ sin 9 }i 

sin - cos <pji 0 

sin Oji cos ^ v sin # 7( sin (p }i cos # 77 y 



(69)(72} 



where dji and ^ are the polar and azimuthal angles of Sphere j relative to Sphere i and are given 
by 



cos Op - (Z f Z,)/fy cos Oji~(Zj- Zj)/rji 



(70a)(73a} 



and 



Pji-timWj-Yd/UCj-Xd] 



This Euler matrix allows the preferred implementation to obtain the &-th sphere's position in the 
primed coordinate system as 



r ki ' = A*r ki . 



mm 



[00261] The current sphere's initial position in the primed coordinate system is given in 
like manner: 



r,,-' = A»r„- 



mem 
mm 



[00262] Now that the three components of the vector r ci ' have been obtained, the 
preferred implementation can find the initial azimuthal position 



#> c , ' = tan" '(jVe, '/*«') 



man 
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If (p c i' is positive, the current sphere will roll toward higher positive values. If it is negative, it 
will roll toward more negative values (higher in magnitude). In the unlikely case that it is zero, a 
random number generator chooses which direction it will roll. 

[00263] The polar angle of the center of the current spheres in the primed coordinate 
system is constant as long as it is in the roll corridor. It is obtained from the law of cosines as 

cos9 c ' = [rji 2 + (a c + a,) 2 - (a c + aj) 2 ][2 r fi (a c + a/)]" 1 m)(7S) 

with 

sinGc' = +[1 - cos 2 G c '] 1/2 . m)(79) 

[00264] There are several possible outcomes as the current sphere descends its roll 
corridor. It can 

• hit the catch net or water level; 

lose contact with Sphere i or Sphere j or both simultaneously; 

• hit the cylinder wall; or 

• hit another pack sphere or several pack spheres simultaneously. 

[00265] The roll corridor ends when the current sphere encounters the first of these 
situations. The preferred implementation then finds the <p c ' at which one of these events happens. 
The event with the smallest magnitude for q> c ' will be the bottom of the roll corridor and will be 
denoted (p c f. 

Hitting the catch net or water leve lH itting the catch net or water level 
[00266] First, the preferred implementation finds the global altitude of the current sphere 
as a function of q> c ': 
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where the inverse Euler matrix is equal to its transposed 4 ~ A* - A' 1 = A T b ecause of its 



orthogonality. The current sphere's global Z component is given by 



Z c = 



Zi + (a c + cc/)(sin8ji sin0 c ' coscp c ' + cosOji cos0 c '). 



Its south pole hits the catch net if Z c = Z net + a c . Solving for the primed azimuthal angle (p c f at 
which this equation may be satisfied, we have 



where the sign of cp c ' is the same as the sign of (p ci \ the angle at the top of the roll corridor. If this 
equation does not have a real solution in (Ojr/2), TO, tc/21, then the current sphere does not reach 
the catch net in this roll corridor. If cp c ' is outside the range \02n), [0, 7t/2"|, then units of 2k are 
either added or subtracted from cp c ' until it does lie within that range. Because neither Gji nor 9 C ' 
can ever be zero or tt, Eq. (79) Eq. 82 is never singular. 

[00267] The preferred implementation also can find which water level the current sphere 
might hit. Water level altitudes between adjacent cylinder walls are generally different as the 
pack is constructed. The only condition under which the current sphere might contact the A:-th 
water level Z w (k) is if and only if the distance P c of the current sphere from the global Z axis lies 
between cylinder wall radii Wk-i and Wk at the point of contact. P c is given by 



cos (p c ' = 



Z ne t + Q c Z x 

a c + a t 



- cos Op cos 9 C 



mm 



sin #, 7 sin 6 C X 



P C = (X C 2 + Y C 2 ) 



{80X83) 



where X c and Y c can be obtained from Eq. (77) Eq. 80 as 



X c =Xi + (a c + ai)(-cos0ji cos^ y/ - sin# c ' cos^ c ' + sin^ sin# c ' sin^ c ' 
+ sinOji cos<pji cosOc) 



(8+a)(84a) 



and 
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Y c = Yi + (cc c + ai)(-cosQji sin^y, sin# c ' cos^ c ' - cos#> 7/ sin# c ' sin^ c ' 

+ smdji sirups cos(9 c ') (S4k)(84b) 



[00268] Using the square of P c rather than Eq. (80), Eq. 83, the water level checked for 
contact is Z w (k)ifW k 2 _ l <P* <Wl . Actually, each annular region between the (fc-l)-th and 
£-th cylinder walls is checked to see if the roll corridor has a possibility of contacting the water 
level there. 

[00269] The expression for the <p c ' at which the south pole of the current sphere touches 
the A>th water level is the same as Eq. (79) Eq. 82 with Z net replaced by Z w {k). Then the end of 
the roll corridor (p c r is either ±x/2 or a contact point with the catch net or a water level as given 
b vEq.(79). Eq.82. 

Losing Contact With One or Both Sphcrcs L osinz contact with one or both 

spheres 

[00270] The preferred implementation also can determine the <p c ' at which the current 
sphere will lose contact with either Sphere / or j or both if it encounters nothing else as it 
descends its roll corridor. 

[00271] Of the two spheres, i and j, the one whose center is higher by the subscript u and 
the one whose center is lower is denoted by 1. Their global positions are (X u , Y u , Z u ) and (X h Yi, 
Zi) and their radii are a u and a u respectively. We always have that Zu > Z/. Consider the 
unprimed coordinate system whose origin is in the center of the lower sphere and whose axes are 
parallel to the global axes. Let 9 C and <p c be the polar and azimuthal angular position of the 
current sphere in the unprimed local coordinate system centered in the lower sphere. 

[00272] As before, the primed coordinate system is the system obtained by tilting the 
unprimed z axis into the center of the upper sphere and then rotating the xy axes of that tilted 
coordinate system about its z' axis until its x' axis is pointing as steeply upward as possible, i.e., 
has a maximum positive projection onto the unprimed z axis. This rotation is what Eq. (69) 
Eq. 72 does if i is replaced by / and j by w. The polar and azimuthal angular positions of the 
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current sphere in the tilted coordinate system are 0c' and (pc' 5 respectively where cos0 c ' is given 
b y Eq. (75) Eq. 78 with the same change in indices mentioned in the previous sentence. 

[00273] The current sphere will lose contact with the upper sphere first. If both lower 
and upper spheres are at the same altitude, it will lose contact with them both simultaneously. 
To find the cpc' position at which the current sphere loses contact with the upper sphere, the 
preferred implementation finds the two longitud e longitudinal lines of the lower sphere in the 
unprimed coordinate system along which the current sphere would just nick the upper sphere if it 
were descending along those longitud e longitudinal lines. The two cp c at which that happens are 
the ends of the roll corridor at which the current sphere would lose contact with the upper sphere 
and roll on the lower alone. Because the current sphere is rolling down only one side of the roll 
corridor, only that side's <p c is of interest. 

[00274] In the unprimed local coordinate system, the current sphere is at 



r c = {&c + sin 6 C cos 6 C + j sin 6 C sin <p c + k cos 9 C ) 



(82X85) 



and the upper sphere is_at 



r u = r u (i sin O u cos 0 U + j sin 6 U sin <p u + k cos 6 U ) 



mm 



[00275] Consider the equation 



r c - r u \ = a c + a u . 



(84)(87) 



[00276] If both sides are squared and written explicitly, we obtain 



[(a c + ai) 2 sin0 c coscp c - x u ] 2 + [(a c + a?) 2 sin0 c sincp c - y u ] 2 + 
[(a c + af) 2 cosGc - z u f = (a c + a u f 



(85aX88a) 



which can be rewritten 
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— — — — — = x u sin 0 C cos ^ c + j/„ sin 0 C sin ^ c + z u cos # c . (85b)(88b} 

2{a c + a,) 



[00277] The components jc M) j/ mj and z w are given in Eq. (83). Eq. 86. If they are 
substituted into Eq. (85b) Eq. 88b and both sides are divided by r Uy we obtain 



+ (®c + a/) 2 - (a c + a u ) 2 . _ . - 

— = sin O u cos g> u sin 6/ c cos <p c + 

sin sin #> u sin 0 C sin + cos 0 M cos 0 C 



where we recognize the l e ft hand left-hand side as cos9 c ', the polar angle in the primed 
coordinate system which is constant as the current sphere descends the roll corridor. Now-Eqr 
{%€ Y Eq. 89 can be written as 

cos9 c ' = sin9 u sin9 c (coscp u cosq> c + sincp u sinq> c ) + cos0 u cos9 c (87a) (90a) 

which through a trigonometric identity can be written 

cos9 c ' = sin9 u sin9 c cos((pc - <p u ) + cos9 u cos0 c . (87b) (90b) 

[00278] We wish to solve for 9 C from this equation. For a given cp c , the 9 C solution has 
three possibilities: 

[00279] 1. There are two 9 C solutions if the current sphere hits the upper 
sphere as it rolls down the cp r longitude longitudinal line of the lower sphere - one solution 
where it touches the upper sphere from above and the other where it touches from below; 

[00280] 2. There is one solution (the two solutions having become 
degenerate) where the current sphere just nicks the upper sphere while rolling along the (p c 
longitud e longitudinal line; and 

[00281] 3. There are no real solutions, meaning the current sphere 
altogether misses the upper sphere while rolling along the cp c longitud e longitudinal line. 



70 



[00282] The second possibility is the one of current interest because that is the point at 
which the current sphere can lose contact with the upper sphere and roll on the longitud e 
longitudinal line of the lower sphere. We will call that point the "departure point." Therefore, 
the preferred implementation first finds the angular position (0 C , (p c ) of the departure point and 
then finds the azimuthal angle (p c ' in the primed coordinate system to which this angular position 
corresponds. 

[00283] To solve for d c , given a particular <p Cy the preferred implementation recasts 
Eq. (87b) is r e cast Eq. 90b into the form 



cos 8 r x 

= sin a sin 0 C + cos a cos 9 C (88)(91) 



^/sin 2 6 U cos 2 {p c - <p u ) + cos 2 9 U 



where 



sin 6 u co$((p c - <p u ) ro/N w ^ x 

sina^ KYc Yu) (89a)(92a} 

^/sin 2 9 U cos 2 (<p c - cp u ) + cos 2 9 U 



and 



cos 9 U 

cosa^-p^ (g9fe)(92b) 
^/sin 2 0 U cos 2 (<p c - <p u ) + cos 2 9 U 



[00284] Applying our trigonometric identity to the right hand right-hand side of-Eqr 
fS& VEq. 91, we obtain 



cos 9 ' 

cos(fl e -g)= c (90)(93) 

>/sin 2 0„ cos 2 - + cos 2 #„ 



which gives the solution 
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9 r = OL + COS 



cos 9 C ' 

^/sin 2 # u cos 2 (#> c - + cos 2 



[00285] The arccosine generally has two distinct solutions: (1) the principal value 
which lies in the interval-fO^h [0, n) and which is denoted by capitalizing the function as Cos" 1 
and (2) the negative of the principal value. We are interested in the case where these two 
solutions are degenerate. That can only be so if the arccosine and its negative are the same 
value, i.e., zero. Therefore 9 C - a at the departure point, and the argument of the arccosine in-B% 
(944- Eq. 94 must be unity: 



cos e; 



■sjsin 2 6 U cos 2 (<p c - <p u ) + cos 2 6 U 



= 1. 



(92X95) 



[00286] This equation is the condition used in the preferred implementation to find the 
unprimed azimuthal position of the departure point. Squaring both sides gives 



cos 2 9 c '= sin 2 0 u cos 2 ((p c - cp u ) + cos 2 0 u 



mm 



which can be <p c to yield 



<Pc = <Pu + COS 



-1 



/ cos^'-cos 2 ^ 
sin : 



(94X97) 



[00287] The square root has a positive sign because the difference between (p c and q> u 
can be no greater in magnitude than 7i/2. Again, the arccosine has two solutions: the principal 
value Cos" 1 and its negative -Cos" 1 . The current sphere would lose contact with the upper sphere 
if it is at azimuthal position cp u ± Cos" 1 . Its true position depends on where it was initially. If its 
initial position (p ci > <p u , then the principal value Cos' 1 is used. If its initial position cp ci < cp u , then 
the negative of the principal value is used -Cos" 1 . In the unlikely case that the initial position 
equals (p U5 then the current sphere is balanced on the upper and lower sphere and a direction is 
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randomly chosen by the random number generator. The preferred implementation using these 
relationships is_to determine where contact is lost. 

[00288] The preferred implementation also determines (0 C? (p c ) of the departure point. It 
then finds the departure point in the primed coordinate system by use of the Euler transformation 
o f Eqs. (69) and (73), Eqs. 72 and 76. 

[00289] The first component of the primed vector is given by 



sin0 c ' cos(p c ' = -cos6 u sin9 c (cos(p u coscpc + sincp u sincp c ) + sin6 u cos0 c (9$) (98) 



where the expression in parentheses can be written cos(cp c - cp u ). Eq. (9 4 ) Eq. 97 can substitute 
for cos((p c - cp u )- Also because 9 C is equal to a o f Eqs. (89), Eqs. 92a and 92b, we have 



. Jsin 2 9 u -sin 2 0 c 
sin 0 C = — 



cos 0 C 



(96a)(99a) 



COS 0 r = 



cos 0 U 
cos 0 r ' 



(96b)(99b) 



[00290] Substituting Eqs. (94) and (96) into (95), Eqs. 97 and 99a/99b into 98, we 
obtain for the end of the roll corridor due to departure 



, cos 0 U sin 0* tan#' 

cos <p c / = = - 

sin 0 U cos 0c 1 tm0 u 



m )(\00) 



[00291] The preferred implementation uses this relationship to obtain 0 C and <p c . Notice 
that the current sphere cannot be in stable contact with two pack spheres unless 9 C ' < 6 U , 
Therefore <pcf always lies between 0 and n/2. Also, remember that cos0 c ' is given by the-4eft 
haftd- left-hand side o f Eq. (86) Eq. 89 and that sin#c' = +[1 - cos 2 9 C '] V2 . 

[00292] The second component of the primed vector is given by 
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sin9 c ' sin(p c ' = (sin(p u coscp c - cos(p u sinq> c ) sin0 c 



WOOD 



and the same set of substitutions leads to 



sin <p C f = 



^/sin 2 0 U - sin 2 6 C 



(99)(102} 



cos # c 'sin O u 



[00293] The <p cf ' defined b y Ego. (97) and (99) Eqs. 100 and 102 is the departure point 
and is the end of the roll corridor. At this position, the current sphere loses contact with the 
upper sphere (or both if their centers are the same altitude). These relationships also are used in 
the preferred implementation to determine 6 C and (p c . 



[00294] The cylinder wall for the current sphere is located at global radius W. The 
preferred implementation can find whether the roll corridor of the current sphere on the two pack 
spheres is such that the current sphere's center can intersect the cylinder wall. If it does, then the 
(p c ' of intersection is found. 

[00295] We need the global position of the current sphere as a function of cp c '. Again, 
the global position of the current sphere is 



where R { is the global position of one of the contact spheres and r c is the current sphere's 
position in the unprimed local coordinate system. (Remember the unprimed coordinate system is 
centered in Sphere i and has axes parallel to the global axes.) The unprimed vector is given in 
terms of the primed vector (that contains cp c ') by the inverse Euler matrix equation 



Hitting The Cylinder Wal lH ittinz the cylinder wall 



R c = R { + r c 



r c = A-'.F c ' 
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where^4"^ - A' 1 is the transpose of the Euler matrix in Eg. (69) and r j - Eg. 72 and rj = (x c ', y c ', z c ') 
is given by 

Xc = (a c + ai) sin9 c ' coscp c ' (102a) (105a) 

y c ' = (a c + a,) sin9 c ' sin <p c ' (102b) (105b) 

z c ' = (a c + a,) cos9 c '. (102c ¥ 105c) 

[00296] The preferred implementation therefore finds the global position-Jk- R c = (X c , 
Y c , Z c ) in terms of the single variable (p c ' describing the descent down the roll corridor. If the 
lateral global position of the current sphere reaches the cylinder wall radius W for a value of (p c ' 
in the interval ((p C j' (p C f') 5 then the current sphere can reach the cylinder wall if it does not 
encounter another object first. In other words, if there is a solution of 



x c 2 + Y c 2 = w 2 m &am 



for 9 C ' in the interval ((p cj ', <p cf ' ), then the current sphere reaches the wall. Written explicitly, 
Eq. (103) Eq. 106 has the form 

A 0 + A] cos(p c ' + A 2 cos 2 (p c ' + Bj sinq> c ' + B 2 sin 2 <p c ' = 0 (4Q4)(!07) 



where 



A 0 = (X 2 + Y 2 - W 2 l{a c + aft + sinO,,- cos9 c ' [2X ( coscp y , + 

2Yi sin(p y , + (a c + a,) sinfy cos9 c '] (105a) (108a) 

A, = -2 sin9 c ' cosO,-/ [X ( cosq),, + Y t sin(p y/ + (a c + a,) sinG/, cosG c '] (±Q5fe)(108b) 
A 2 = (a c + a,) cos 2 6/, sin 2 9 c ' (WSe)(i08c) 

B, = 2 sin9 c ' [AJ sincpy,- - Y t cosq),-,] (4QM )(108d) 
£ 2 = (a c + a,) sin 2 G c '. (W§a)£108e) 
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[00297] Eg. (10 4 ) Eg. 107 can be converted into a guartic eguation in sin(p c ' by changing 
the cos 2 (p c ' term into 1 - sin 2 (p c ' and the coscp c ' term into +[1 - sin 2 <p c '] 1/2 . (The plus sign only is 
used only because (p c f cannot exceed ±rc/2.) We then have 

(A 0 + A 2 ) + Bj sin(p c ' + (B 2 - A 2 ) sinV - -A } [\ - sin 2 (p c '] 1/2 . <4O6¥109) 

[00298] Sguaring both sides gives 

(A 2 - B 2 ) 2 sinV - 2Bj(A 2 - B 2 ) sinV + {A 2 - 2AqA 2 - 2A 2 + Bj 2 + 2AqB 2 + 

2^2^) sinV + 2(^ 0 + ^2)5/ sin^ c ' + (A 0 -Aj+ A 2 )(A 0 + + ^2) = 0 (407 X110) 

[00299] If this eguation has any real roots in the interval ((p c /, <p C f'), then the smallest in 
magnitude of these that satisfies the original e guation (10 4 ) Eg. 107 is the point at which the roll 
corridor intersects the cylinder wall. The preferred implementation can check the roots in the 
original eguation because sguaring can give extraneous roots. 

Hitting One or More Spheres Simultaneousl yH itting one or more spheres 
simultaneously 

[00300] The preferred implementation has now defined a roll corridor in the primed 
coordinate system beginning at <p ci ' and ending at q> c /. The end of the roll corridor is where the 
current sphere either lost contact with one or both of its contact spheres, when its south pole 
encountered either a catch net or a water level, or when its center encountered its cylinder wall. 
The preferred implementation can check all spheres that can intersect this roll corridor to find the 
sphere(s) encountered first as the current sphere descends along its roll corridor, i.e., as <p C ' moves 
from (pd to (p c /. 

[00301] The global position of the current sphere as a function of (p c ' is given b y Egs. 
(77, 78, and 81). Egs. 80, 81, and 84. If some other pack spheres, say the k-th sphere is to 
intersect the roll corridor, then the following eguation must have a solution in the range (<p ci \ 

<PcfY 
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R c (<p c ')-R k =a c + a k . (W8)0II) 
[00302] Squaring this equation gives 

Rl - 2R C • R k + Rj = (a c + a k f am(U2) 
which becomes the following equation in <p c ': 

A cos<p c ' + B s\nq> c ' = C ^4^ (113) 

where 

A = -X ik cos <pji cos dp sin 6 C ' - Y ik cos dp sin <p Jt sin 6 C ' + Z ik sin d c ' sin 0jj(444X114) 

5 = X ik sin pj f sin B c ' - Y ik cos ^ sin dc' (443 4(115) 

C = -A!}* cos cos 0 C ' sin dp - Y ik cos d c ' sin p,,- sin - Z ik cos 0 C ' cos dp 

(2a c + a, + a k )(a k - a.) - Rl 

+ — „ * A x ,; (445)016} 

2(fl c + a,) 

[00303] The global coordinates of Sphere k relative to Sphere i are given by 

Xik = Xi - X h Y ik = Yi - Y k , Zik = Zi - Z k , 044aXil7a} 

and the square of the distance between the two sphere centers is 

Rik = X ik 2 + Y ik 2 + Z ik 2 . (444fe)(117b) 

[00304] The solution procedure o f Eq. (110) Ea. 113 is similar to that given in Section 
VII.D.7. The preferred implementation finds <p c ' to be 
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<p c ' = a± cos" 1 , C r m5e ¥118a) 



where 



a = tan" 1 -. (444b)(U8b} 
^4 



[00305] If any of the <p c r solutions are outside the range [0,270, [0, 2tz]^ units of 2n are 
added or subtracted until they are inside that range. Then if any solutions lie within-^-^/) 
( (pr[, (pr/ ) the smallest in magnitude of these is the position of first contact between the current 
sphere and Sphere k. If there is no solution, then Sphere k intersects no part of the roll corridor. 

[00306] The preferred implementation has now performed an exact solution to find 
whether and where Sphere k intersects the current sphere's roll corridor. In most packs, the vast 
majority of the pack spheres are out of range of the roll corridor and so do not need the detailed 
calculations o fEqs. (115). Eqs. 118a and 118b. There are a few much simpler checks that can 
be performed by the preferred implementation to immediately determine whether a pack sphere 
is out of range and so is not a contact candidate. 

[00307] First, the roll corridor's highest global altitude is at <p ci ' which, from Eq. (78) 
Eq. 81 is found to be 

Z c ,max = Zi + (a c + a,-)(sin0ji sin0 c ' coscp ci ' + cos0jj cos0 c '). (H6) (ll9) 

Therefore, the k-th pack sphere whose center lies a distance a c + ak above Z C)tnax cannot overlap 
the roll corridor. Likewise, the lowest global altitude is at q> c f and is given by 

Z c ,min - Zi + (a c + cc/XsinOjj sin0 c ' cos(p cf ' + cosfy cos0 c r )- (4434 (120) 

[00308] The A>th pack sphere whose center lies a distance a c + ak below Z Cymin also 

cannot overlap the roll corridor. Therefore, the preferred implementation can immediately 

exclude from consideration all spheres whose global altitudes Z k satisfy the following: 
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k-th sphere does not overlap if 

Zk > Zc^ox^ a c + at 



or 



Zk < Z Cy min ~ 0. c - <Xk. 



[00309] Roll corridor extrema in the X and Y directions do not necessarily lie at the 
beginning and end of the roll corridor. Bqs. (81) Eqs. 84a and 84b g ive the global X and Y 
coordinates of the current sphere as a function of <p c \ If we differentiate X c (<p c ') with respect to 
<Pc and set the derivative equal to zero, we have an equation whose solution gives us the extrema 
in X. If either of those extrema lie 'm4pJ^<f± ((pn, (prf l then that position is either the minimum 
or maximum X value the current sphere assumes as it descends the roll corridor. The same is 
true for the derivative of Y(<p c ') with respect to <p c '. The derivatives followed by the extrema 
solutions are 



dX c 



- = (a c + a,)(cos Op cos (p }i sin 9 C 1 sin <p c 1 + sin <p yi sin 9 C 1 cos <p c ') = 0 (119 X 122) 



d(p c 



tan^' = - 



tan (pji 

COS0.; 



for the two X extrema (the above equation has two solutions for <p cx '), and 



dY c 



- = (a c + a,)(cos e n sin <p M sin 0 C 1 sin <p c ' - cos <p }i sin 9 C 1 cos <p c f ) = 0 (4344 (124) 



d<p c 



tan^ f = - 



cot <Pji 
cos 9 n 



fm) (125) 
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for the two Y extrema. 

[00310] Finally then, the X and Y extrema of the roll corridor in the global coordinates 

are 

X c , min = MN[X c (<p ci % X c {(p c f), X c (<p cx ')} fl33a ¥126a) 

X c , max = MAX[X c (<p ci '), XM> c f), XApac 1 )] (t33fe)026b) 

with the X c (<p cx ') omitted if neither yd lies 'vaAtyJ^pefh - ((pn\ Similarly, 

Y c>min = MESf[y c (9) c/ '), Y c (<p c /), Y c (<p cy ')] (i34a)027al 

Y c , max = MAX[Y c (<p ci '), Y c (y c /), Y c (<p cy ')], mm iUlb) 

again with the Y c (<p cy ') omitted if neither (p cy ' lies in-tyJ^p^ (cpr/. <pr/ ). 
[00311] Then the k-th sphere cannot overlap the roll corridor if 

X k < X cMn -a c -a k <435a ¥128a) 



or 



X k > X Ctmax + a c + a k (±2$fe)Q28b} 



or 



Y k < Y c/nin - a c a k <435e 4(128c) 



or 



Y k > Y c , max + a c + a k . <43M)02Ml 
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[00312] Another discriminator that can be used by the preferred implementation to 
immediately exclude a large number of pack spheres from detailed examination is to consider 
whether the roll corridor is far enough from the central Z axis that it cannot overlap all the 
spheres within a given inner cylinder. To find this, the preferred implementation uses the closest 
radial approach P Ctmin of the roll corridor to the central axis. The square of the radial distance 
[P c ((Pc)] 2 is given by 

[Pc(<Pc')? = [X c {<Pc')i + [ W)] 2 (136)029} 

where we obtain the explicit expression by substituting Eqs. (81) into Eq, (126): Eqs. 84a and 
84b into Eq. 129: 

[P c {<Pc)\ = [Xi + (a c + ai)(-cos^7 cos^/7 sin# c ' cos<p c ' + sin^y/ sin# c ' sin#> c ' 

+ smdji costpji cos # c ')] 2 
+ [Yi + (a c + ai)(-cos6ji sin^ 7 sin# c ' cos<p c f - cos^y,- sin# c ' sin#> c ' 

+ sindjt sintpji cos# c ')] 2 . f±2£ K130) 

[00313] Setting the derivative of [P c (<pc)f with respect to <p c ' equal to zero yields the 
equation whose solution (p cp f gives the closest approach of the roll corridor to the global Z axis: 

[Xi + (a c + ai)(-cosdji cos^ sin#c' cos#> c ' + sin^-/ sin# c ' sm(p c ' + 

sin^v cos^y,- cos# c ')](cos#y, cos<pji sin8 c ' sm<p c r + sin^y,- sin# c ' cosp c ') + 

[7, + {a c + aj)(-cos6ji sin^y, sin# c ' cos<p c ' - cos^y,- sin#c' sin^ c ' + 

sinffji sin^y; cosfl c ')](cos/9y; sinpy/ sinflc' smg) c r - cospy, sin# c ' cos#? c ') = 0.R284 O31) 

[00314] This equation is of the form 

A Q +A } cosp c ' + A 2 cos V + Bi sinp c ' + B 2 sinV = 0 f!39 4(132) 

where in this case, the coefficients are given by 
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Aq = X 2 + Y 2 + 2(a c + fl/) cos# c '(^7 sin^/ cos^y + 7/ sin#,y sin^y) + 

(a c + a,) 2 cos 2 # c f sin 2 0 y7 f!30a) (133a) 



^4/ = -2(a c + a t ) cosOji sinft/ [Xi cos^y + 7/ sin^y + 

(a c + a,) cos# c ' sin^y] fl3Qb ¥133b) 

^ = (a c + af cos% sin 2 0 c ' ft3Qe ¥133c) 

,87 = 2(a c + a,-) sin# c ' {smdji [X t + (a c + ai) cos# c ' sin0 y y cos^y] 

- c6s6 c '[Yi + (a c + a/) cos# c ' sin0 y/ sin^y] } (±30d)(133d) 

B 2 = (a c + a/) 2 sin 2 # c ' (i50e)£133e) 

[00315] Eg. (129) Eg, 132 can be converted into a guartic eguation by changing the 
cosines into sines in the same fashion as led up to Eg. (107): Eg. 110: 

(A 2 - Bif sinV - 25/042 - B 2 ) sinV + (Aj 2 - 2AoA 2 - 2A 2 2 + 
Bj 2 + IA0B2 + 2^ 2 5 2 ) sinV + 2(A 0 + A 2 )B } sm(p c ' + 

-A } + + ^7 + A 2 ) = 0. fm4(134) 

[00316] There are four sin g> c f solutions to this quartic equation and for a given sin <p t r 
there are two <p c ' solutions. Complex solutions and real solutions lying outside \ 1,1] [1, 11 are 
immediately discarded in the preferred implementation. Furthermore, we are only interested in 
those solutions with the same sign as sin#> c / and that lie in-(y^-y^9- {(f> c j, (prf ) because that is the 
direction the current sphere will roll. The preferred implementation checks any remaining 
solution (p C p by ensuring it satisfies the original equation (129). Eq. 132. If there is no solution 
in-^J^f^ (y?w\ (prf ) then one or the other of the end points is the closest approach to the global 
Z axis and both must be substituted into Eq. (127) Eg. 130 to find the smaller P c . 
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[00317] Once the P C jnw of closest approach is found, the preferred implementation 
checks to see if there are any mode's cylinder radii Wk that lie completely inside P c ,mm by more 
than the sum of the mode sphere's radius plus the current sphere's radius. If 

P c >W k + a!+a c fl32 ¥135) 

then none of the Mode k particles can touch the current sphere as it descends its current roll 
corridor and so may immediately be excluded as candidates for contact. 

[00318] All spheres whose north poles have sunk beneath the catch net or water level are 
not candidates for contact with the current sphere. The preferred implementation continually or 
periodically monitors the lowest index number of the sphere whose north pole is still above the 
catch net so that w e wast e no time is wasted looking at spheres with lower index numbers as 
candidates for contact. 

Touching One Sphere and The Cylinder Wal l Touching one sphere and the 

cylinder wall 

[00319] Remember that in the preferred implementation it is the center of the current 
sphere constrained by the cylinder wall, not its edge. When in stable contact with both a pack 
sphere and the cylinder wall, the roll corridor lies on the wall and descends toward the equatorial 
plane of the contacted pack sphere. If the contacted sphere (denote it Sphere i) is as large as the 
cylinder, then the current sphere moves toward the lowest point on a roll corridor that spans the 
entire 2n radians of the cylinder. 

[00320] The roll corridor will be denoted by the global azimuthal position 0 C . Its initial 
position is given by 

& c0 - tan" 1 (YJXc). (133)036) 

[00321] The preferred implementation uses the range r0,27i) TO, 2k] rather than ( n,n] 
r-7c, ri\ for the range of the global azimuthal position. The global azimuthal position of the 
Sphere / is 
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&i = tan 1 (Yi/Xi) 



(144)037} 



which also lies in [0,2ft). [0, 2k], If & c o > 0i the current sphere rolls toward higher <P If 
0 c o < 0h it rolls toward lower 0. In the unlikely event that 0 c o = 0» the random number 
generator chooses a roll direction. 

[00322] If it encounters nothing else, the current sphere will roll until it reaches: 

[00323] 1. the equatorial plane of Sphere / where it will fall away from 

Sphere /; or 

[00324] 2. the lowest azimuthal point in its cylinder permitted by Sphere i. 
[00325] The azimuthal position at which it will reach the equatorial plane of Sphere / is 
found by the law of cosines as 

0 cf = 0 i ±A0 c msxns) 



where 



Ac2> c = cos _1 



' W?+P>-(a e + ai )» 

v 2*y>, J 



[00326] In this expression, W c is the radius of the cylinder wall, P, = ^jxj 1 + Yf is the 
radial position of Sphere / from the global Z axis, and a c and a t are the radii of the current sphere 
and Sphere i, respectively. The plus sign is used in Eq. (135) Eq. 138 if the current sphere is 
rolling toward increasing 0, the minus sign is used if the current sphere is rolling toward 
decreasing 0. If the absolute value of the quantity in the parentheses o f Eq. (136) Eq. 139 is 
greater than unity, then there is no solution and the current sphere simply rolls around the 
perimeter of the cylinder until it reaches the lowest point it can. This lowest point is opposite the 
azimuthal position of Sphere i. Hence if 
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w c 2 +p?-(a c + ai y 



2W c P t 



>1 



then 



A@=±n 



[00327] In the unlikely event that P t = 0, i.e., Sphere i is centered on the global Z axis, 
then there is no downward direction. The roll corridor is horizontal and the current sphere is 
placed where it is. 

[00328] In this analysis, if the addition or subtraction of A<P C causes @ C f to lie outside 
[0,2ft), [0, 2n], then units of 2n are added or subtracted until 0 C / again lies in \0,2n). TO, 2%]. 

[00329] The preferred implementation also determines what placed pack spheres might 
intersect the roll corridor in the range-f^o^^- ( 0 ^, 0^ . 

[00330] The current sphere will impact a placed pack sphere, say the >th sphere, if there 
is a solution to 



R c (O e ) - Rj =a c + aj 



anywhere in the range-4^w*4yh- (flw ,. 0 r J i. The Cartesian components of the vector 
R c (O c ) = X c , Y c , Z c are given by 



X c (0 c ) = W e COS0 c 



(M0)£143) 



Y c (0 c ) = W c sm0 c 



(441 4(144) 



[00331] The Z c component can be found by inserting the X c and Y c component into the 
condition 
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\R e -R,=a e + a, . (442 4(145) 

[00332] Squaring both sides, we obtain 

W 2 + Zc + R 2 - 2 W c (Xi cos@ c + Yi sintf> c ) - 2Z,Z C = (a c + a,) 2 (W3)£146) 

which is quadratic in Z c . One of the solutions is when the current sphere is above Sphere i and 
the other when it is below. In the preferred implementation, we are only interested in the 
solution above Sphere i, which is 



Z C (0 C ) = Zi + [Z- + (a c + a,) 2 - W 2 - R t 2 + 2W c (Xi cos0 c + 7, sin0 c )] 1/2 .(M4)Q47) 



[00333] Now, insertin g Eqs (1 4 0. 1 4 1. Ml), in Eq. (139). Eqs. 143. 144. 147 in 
Eq. 142, we obtain an equation to solve for the <P C at which the current sphere contacts the 7-th 
pack sphere: 

A 0 + A, cos&c + B, sm& c + (A 2 cos0 c + B 2 sin# c ) 2 = 0 (M5)£148} 
where the coefficients are given by 



A 0 =G 2 - 4(Z, - Z y ) 2 [Z, 2 - W c 2 - R 2 + (a c + a,) 2 ] <446a ¥149a) 



Ai = 4G W c {Xi - Xj) - 8(Z, - Zj) 2 W c Xi &46b )(149b) 



A 2 = 2 W c (Xi - Xj) fM6e ¥149c) 

B, = AGW c (Yi - Yj) - 8(Z, - Zj) 2 W c Y t (446d ¥149d) 

B 2 = 2 W c {Yj - Yj) &46e ¥149e) 
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and G is given by 



G = (at - aj)(2a c + a t + aj) - R 2 + /?/ + 2Z/ 2 - IZiZj . (443 X150) 

[00334] Eq. (1 4 5) Eq. 148 can be converted into a quartic polynomial in cos0 c : 

{A 2 + #2 2 ) 2 cos 4 0 c + 2(^7^/ + 2A 2 BjB 2 - ^y5 2 2 ) cos 3 <Z> c + 
(Aj 2 + + 5 7 2 - 2AoB 2 2 - - 2B 2 4 ) cos 2 & c + 

2(/M, - 2A 2 BiB 2 +AiB 2 ) cos0 c + 

- 5 y + + 5 y + 5/) = 0 (448¥1S1) 

whose solutions are candidates for contact of the y-th pack sphere. The preferred implementation 
uses this expression as part of its determination. There are four solutions of the quartic equation 
and each cos^ c solution has two 0 C solutions — eight in all. Of the four quartic solutions, the 
preferred implementation can eliminate complex solutions or solutions whose magnitude is 
greater than unity. Of the remaining solutions, if there are any that satisfy the original e quation 
<44SV Eq. 148 and that lie in the range of the roll comdor-f^a^^r ( 0^ they are contact 
points. Of the contact point(s), the value of 0 C that lies closest to 0 c o is the point at which the 
current sphere touches the y'-th sphere. We will denote that point by & cj . If there are no solutions 
that satisfy the above conditions, the current sphere cannot contact the j-th sphere. 

[00335] As the preferred implementation checks the candidate spheres in the pack for 
contact with the current sphere, that sphere with 0 cJ closest to 0 c o, if any, is the next sphere 
contacted. In the unlikely case that two or more pack spheres have the same 0 CJ , then it checks 
the stability of each trio of contacts to see if the current sphere is to be placed there. If it is not 
stable, then the stability routine determines which sphere(s) the current sphere will roll away 
from. 

[00336] As an efficiency measure, the preferred implementation can eliminate from 
consideration any modes or categories that are sufficiently far inside the cylinder wall of the 
current sphere that they cannot touch it. All mode k spheres are out of range if 
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W k + a k <W e -a e . 



(±49)(152) 



[00337] Any other sphere j is a candidate for being in range only if 



W c -a c - aj < Pj <W c + a c + ctj 



(440)053} 



where Pj is the global radial distance of Sphere j from the Z axis. 

THE SUBROUTINE UPDAT E THE SUBROUTINE UPDATE 

[00338] The UPDATE subroutine of the preferred implementation updates all 
parameters and arrays that change as a result of the placement of the current particle. These 
changes include the following. 

1. The Pool Switch 

[00339] If the pool switch is true, then the placement of the current particle means the 
bin from which it was chosen should be diminished by one. This adjustment is performed in 
Subroutine UPDATE. 

2. The Catch Net 

[00340] The catch net belonging to the mode, submode or category current sphere also 
generally ratchets upward. The prescription for determining how much it ascends will now be 
addressed. 

[00341] A sense of pack surface must be established. With a variety of particles being 
dropped into a variety of cylinders, the pack surface must be defined. This can be done in a 
number of ways, as noted above. For purposes of the preferred implementation, let us consider 
the altitude of the south poles of the last m placed spheres of any and all modes. As explained 
above, the average of those altitudes is defined for present purposes as the pack "surface": 



'surface 



1 m 

= — I (Zj - Qj) 

m >i 



vn.E.2fn ri54) 
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where m is the number of most recently placed pack sphere equatorial cross-sections that would 
just cover the outermost cylinder n with cross-sectional area nW n 2 . The value for m is found 
from the inequality 

m-\ m 

Z af <W?<Z aj . VII.E.2(2) 055) 

y=i 7=1 

[00342] If the size of the last placed set of spheres is very large so that m is very small, 
then we have a minimum required value for m which can be user-defined but for which the 
default value in this implementation is m min = 12. Only when the pack is just starting so that 12 
particles have not yet been dropped is m smaller for this implementation. Pr e plac e d Preplaced 
particles do not count in assessing the pack "surface." 

[00343] The A>th mode catch net is then defined as 

Znet(k) = Zsurface ' JClk VIIE.2(3) (T 56) 

where y is the number of k-th mode radii beneath the surface that the k-th mode catch net lies. 
Here it is a user-defined number which has a default of y = 2. 

[00344] Subroutine UPDATE readjusts all the above parameters so that the catch net for 
every mode has a proper value after the current sphere is placed. 

[00345] When the pack is still small enough that an insufficient number of particles have 
been placed to satisfactorily define a surface, the catch net may be ratcheted up a little after the 
placement of each particle so that the particles on the bottom of the cylinders do not all have the 
same altitude. If the switch ROUGHNET is true, then the current sphere's catch net moves 
upward by a small random amount 

ra min IN cs VII.E.2( 4 ) (157) 

where r is a uniformly generated random number in KU1 TO, 11 and N cs is the total number of 
sphere cross sections required to cover all the cylinders: 
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N CS =T. 



ft 




158) 



where n is the number of modes, Wk is the cylinder wall of the £-th mode, and Wq = 0. If 
ROUGHNET is false, then the catch net stays where it is until enough particles have been 
defined to make a surface and until the surface is more than Z sur face - y«* above the starting place. 

[00346] The starting place for the &-th mode catch net is Z in i t - a*. A A>th sphere radius is 
subtracted from the overall pack starting altitude Z init because the catch net stops south poles and 
Zinit is the starting place for sphere centers. The default for Z, m/ is zero in this implementation. 

The Water Lcvels The Water Levels 

[00347] The water level also is modified when the current sphere is placed according to 
the preferred implementation. Recall that the water level only exists beyond the innermost 
cylinder. Its altitude tracks the pack surface and is at a user-specified distance P in terms of the 
smallest sphere radius a m i n beneath the surface defined b y Eqs. VII.E.2(T) and (2). Eqs. 154 and 
155. The default for P is unity. Hence 



If P is negative, the water level lies above the pack surface. 

[00348] The smallest cylinder in this implementation has no water level, only a catch 
net. It does not need a water level because water levels are designed to simulate the surface of 
the innermost cylinder and the innermost cylinder has its own particle surface potentially 
containing all the particles of the pack. 

Current Sphere History Paramctcrs Current Sphere History Parameters 
[00349] The preferred implementation uses are parameters that monitor the number of 
times the current sphere tunneled, the number of distinct roll corridors it underwent (not counting 
free-fall drops), and the number of free-fall drops it underwent. All these are reset to zero in this 
implementation. 



7 = 




am 
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[00350] The arrays indicating which pack spheres the current sphere is touching and 
which spheres it just rolled away from are also zeroed out along with the total number of such 
spheres. 

[00351] The array keeping track of the lowest index number above catch net for each 
mode is adjusted if necessary. 

[00352] The array monitoring the index numbers of all pack spheres with exposed north 
poles also is updated if the current sphere covered one or more exposed north poles. 

[00353] In addition, the pack sphere position arrays have the current sphere's (X,Y,Z) 
(X % y, Z) p osition added to them along with the array specifying its mode number. The counter 
monitoring the total number of spheres placed so far in the pack is also incremented by one. 

THE SUBROUTINE OUTSTA T THE SUBROUTINE OUTSTAT 

[00354] There are a number of important pack statistics and features that the user may 
need or want to calculate. In this implementation, subroutine OUTSTAT calculates the output 
statistics that the user specifies. These output statistics include: 
[00355] A. The volume fraction of the pack. 

[00356] B. The radial distribution functions gi/r)dr which is the mean 
number of Mode j particle centers that lie between r and r + dr from a Mode i particle center. 

[00357] C. Coordination numbers Qy which represent the mean number of 
Mode j particles touching a Mode i particle. 

[00358] D. Detailed pack statistics for each particle mode such as that 
mode's minimum and maximum positions, number of its particles in the pack, radial and axial 
density functions for each requested mode, as well as other information. Each of these output 
data will now be addressed. 

Finding the Volume Fraction Finding the Volume Fraction 

[00359] To determine the packing fraction, the program finds the total volume of the 
innermost cylinder occupied by spheres. Recall that the cylinder walls are boundaries to the 
centers of the spheres, not their edges. Hence, many spheres will be overlapping the cylinder 
walls. The preferred implementation has a general expression for determining how much 



91 



volume of a sphere of radius a s located at position (x s , y s , z s ) lies inside a cylinder which is 
concentric to the z axis with radius a c , bottom cylinder altitude z b and top cylinder altitude z t . 
The sphere volume inside the cylinder will be denoted V{ n . Here the sphere and cylinder 
geometric parameters have no restrictions except: 

a s >0,a c > 0, and z, > z b . VIII.A.f 1 X 160) 

[00360] The radial distance of the sphere center from the z axis is defined here by 

Ps = [xs 2 + y s 2 ] m . VIIIA.(2¥ 16n 

[00361] There are four possible cases for this overlap problem. In Case 1, the sphere is 
completely outside the cylinder so V in = 0. In Case 2, the sphere is completely inside the 
cylinder so V in = (4/3)rac s 3 . In Case 3, the cylinder is completely inside the sphere so V in = na c 2 {z t 
- zt). hi Case 4, the sphere and cylinder intersect each other and the expression for V in is more 
complex. In the preferred implementation, processing is performed to identify and address these 
cases using the following approach. 

Case 1: Sphere Completely Outside Cylindc r Cdise 1: Sphere Completely Outside 
Cylinder 

[00362] The sphere is completely outside the cylinder 

V in = 0 VIII.A.Un (162) 
only if one of the following conditions are met: 

If the sphere's north pole is lower than the cylinder bottom ... 

z s + a s < z b Vm.A.l(la) (163a) 
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Or if the sphere's south pole is above the cylinder top 

z s -a s > z t 



)Q63b) 



Or if the sphere is more than its radius away from the cylinder side 

Ps-a s > a c 



)£163c} 



Or if the sphere center is above the cylinder top and is more than its radius away 
from the upper cylinder edge. . . 

(163d) 



(ps - a c f + (z s - z t ) 2 > a s 2 and z s > z t 



Or if the sphere center is below the cylinder bottom and is more than its radius 
away from the lower cylinder edge . . . 

Q63e) 



(p s - «c) 2 + (zb - Zs) 2 > a s 2 and z s < z b 



Case 2: Sphere Completely Inside Cylindc r Case 2: Sphere Completely Inside Cylinder 
[00363] The sphere is completely inside the cylinder 



V in = (4/3)na s 



)064) 



only if all three of the following conditions are simultaneously met: 



If the sphere's north pole is no higher than the cylinder top . 

z s + a 5 < z t 



)(165a) 



And if the sphere's south pole is no lower than the cylinder bottom . . . 

z s -a s > z b Vm.A.2(lb ¥ 165b) 



And if the sphere's equator nowhere exceeds the cylinder radius 
p s + a 5 <a c . 



fl65c) 
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Case 3: Cylinder Completely Inside Sphere Case 3: Cylinder Completely Inside Sphere 
[00364] The cylinder is completely inside the sphere 



V in = na c 2 (z t - z b ) 



Q66) 



if both the following two conditions are met: 



If the upper edge of the cylinder is inside the sphere . . . 

(p s + a c ) 2 + {z t -z s f<a s 2 




Q 67a) 



And if the lower edge of the cylinder is inside the sphere . . . 
(p s + a c f + (z s • z b f < a s 2 . 




067b) 



Case 4: Cylinder and Sphere Intcrsec t Case 4: Cylinder and Sphere Intersect 

[00365] There are several subsets of Case 4. To efficiently organize them, we define 

two circles: the sphere circle and the cylinder circle. The sphere circle is the perimeter of the 

sphere's cross section at a given z altitude. Its radius is 



and its radial position p s is given b y Eq. VIILA(2). Eq. 161. Likewise the cylinder circle is the 
perimeter of the cylinder's cross section which has constant radius a c at any z altitude between z b 
and z t and is centered on the z axis. 

[00366] To describe the intersection of sphere and cylinder, four subcases for the 
relationship of the sphere circle and the cylinder circle at a particular altitude z are defined: 



r s = W - (z - z s f] m , z s -a s <z<z s + a s 




[00367] 



Subcase 4.1: The two circles do not overlap or one or both are 
nonexistent. 

Subcase 4.2: The sphere circle is completely enclosed in the 
cylinder circle. 



[00368] 



[00369] Subcase 4.3: The cylinder circle is completely enclosed in the 

sphere circle. 
[00370] Subcase 4.4: The two circles intersect. 
[00371] To find Vi„, the preferred implementation integrates in z the overlap areas of the 
two circles A in from the altitude where one subcase holds to the altitude where another subcase 
holds. 

For Subcase 4.1, A in = 0. 



For Subcase 4.2, A in - na s 



Vin.A. 4 (2b ¥ 169b) 



For Subcase 4.3, A in = na c 



For Subcase 4.4, A in = a c 2 (a - cosct sina) + r s 2 (fl - cos/? sin/?) 



ri69c) 
f!69d) 



where r s is given b y Eq. VIII.A. 4 (1) Eq. 168 and 



cos a = 



VIII.A.1(3a ¥ l70a) 



sin a = 



yj[r? -(a c -Ps) 2 ][(a c + p s y-rn 
2p s a c 



)(]70b} 



cos P - 



Pi + rj - a 2 c 



2p s r s 



f!71a) 



sin P 



= yJ[r?-(a c -Ps) 2 ][(ac + Ps) 2 -r?] 
2p s r s 



ri71b) 
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[00372] The implementation now constructs a series of altitudes at which the cylinder 
circle and the sphere circle change from one subcase to another in Case 4. 

[00373] A sphere and a cylinder side can intersect each other at either 0, 2, or 4 altitudes. 
Pairs of altitudes may be degenerate (two solutions having the same altitude) if the sphere just 
nicks the cylinder wall. The preferred implementation determines these altitudes using the 
following equation: 

[r s (z)f - [a c ± p s f VIIIA. 4 (5a Y 172a) 

or 

a s 2 - (z - z s f = [a c ± p s f VIII.A. 4 (5b) (172b) 
which, using Eq. (9), Eq. 176, has the four solutions 

z = z,±^ 2 -(a c ±p s ) 2 . VIII.A. 1(6X 173) 

[00374] The + sign in the first ± o f Eq. VIII.A.K6) Eq. 173 denotes the upper or 
northern hemisphere intersection altitude; the - sign denotes the lower or southern hemisphere 
intersection altitude. The + sign in the second ± o f Eq. VIII.A. 4 (6) Eq. 173 denotes intersection 
with the far cylinder wall, on the opposite side of the z axis from the sphere's side; the - sign 
denotes intersection with the near cylinder wall, on the same side of the z axis as the sphere. 

[00375] Armed with these four solutions, the implementation can effectively address 
each of the four subcases of Case 4. In this illustrative implementation, it separately addresses 
whether the sphere circle lies in the southern hemisphere of the sphere or whether it is the 
equatorial plane or in the northern hemisphere. The four subcases are considered separately for 
the two hemispheres making eight situations in all. The hemispheres are considered separately 
because if the sphere circle is in the southern hemisphere, it is increasing in size as z increases. 
Hence the altitude at which we will change subcases can be determined. If the sphere circle is 
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the equatorial plane or in the northern hemisphere, it is decreasing in size as z increases and 
again the altitude at which subcases will change can be determined. 

Finding Upper and Lower Boundaries of Subcases for Case ^ Finding Upper and Lower 
Boundaries of Subcases for Case 4 

[00376] The logic used in the preferred implementation to find the altitudes at which 
subcases change will now be discribed. Let z\ ow be the lower altitude of the subcase currently 
being examined and z^gh be its upper altitude. The function r s (z) is given in Eq. VIII.A. 4 (1). 
Eq. 168. 

[00377] Initially, let z iow = z b , the bottom of the cylinder. After each DO WHILE cycle, 
z iow is set to the previous cycle's z high and a new z high is sought. To illustrate using pseudocode 
expression: 

Vi n = 0 ! Initialize V in to zero. 

Zjow ~~ Z D ! Initialize zi ow to the bottom of the cylinder. 

keepGoing = TRUE ! Do Loop flag. 

DO WHILE (keepGoing) 

IF (z s > ziow) THEN ! Southern hemisphere subcases 

IF (p s - r s (z hw ) > a c ) THEN ! Subcase 4. 1 

IF (zhigh > z t ) keepGoing = FALSE 

(There is no addition to V in in this subcase.) 
ELSE IF Go, + r s {z hw ) < a c ) THEN ! Subcase 4.2 

z high = MINtz, - ^Ja* -{a c - p s ) 2 , z„ z s + a 5 ] 

Vin = Vin + K{Zhigh ~ Zi ow )[a^ - Z s 2 + Z s {Zhigh + Zfow) - 

2 2 
(Z high + ZhighZlow + Z fow)/3] 

IF (zhigh = z t OR Zhigh = z s + a s ) keepGoing = FALSE 
ELSE IF (r s (z, ow ) - p s > a c ) THEN ! Subcase 4.3 

z high = MIN[z, + y/af - (a c + p s f , z t ] 

Vin = V in + Tide (Zhigh - Ziow) 

IF (zhigh ~ z t ) keepGoing = FALSE 
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ELSE IF ((r s (zi ow ) > p s AND r s (zi ow ) -p s <a c AND a c < r s {z hw ) + p s ) 
OR (p s > r s {z, ow ) AND p s - r s (z, ow ) < a c AND a c <p s + r s (z, ow ))) 

THEN ! Subcase 4.4 

IF {a s >a c + ps) THEN 

Zhigh = MW[z s - J a* ~(a c + p s ) 2 , z t ] 

z high 

V in =V in + J dzA in = V in + AVi VIII.A. 4 (7¥ 174) 

Z low 

where AVi is given in the next section below. 
IF (zhigh = z t ) keepGoing = FALSE 

ELSE 

z h i gh = MN[z s + y/af - (a c - p s f , zj 

z high 

V in = V in + J dzA in = V in + AV 2 VIII.A. 1(8 X 1751 

z low 

where AV2 is given in the next section below. 
END IF ! End two Subcase 4.4 situations 

END IF ! End Subcase 4.4 

ELSE ! Equatorial plane & northern hemisphere subcases 

IF (p s - r s (z low ) > a c ) THEN ! Subcase 4. 1 

keepGoing = FALSE 

(Because the sphere circle is shrinking away from the cylinder 
circle, there is no chance for additional overlap.) 

ELSE IF (p 5 + r s (z hw ) < a c ) THEN ! Subcase 4.2 

Zhi g h = MIN[z,, z s + a s ] 

Vin = V in + n(zhi g h - z iow ) [a s 2 - Z? + Z s {z h igh + Zlow) - 

2 2 
(Z high + ZhighZlow + Z / 0 w)/3] 

keepGoing = FALSE 
ELSE IF (r s {z hw ) - p s > a c ) THEN ! Subcase 4.3 

z U g h = MIN[z s + yjaf - (a c + p s ) 2 , z t ] 

Vin = V tn + Tta C 2 (Zhigh - Zhw) 
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IF (zhigh = z t ) keepGoing = FALSE 
ELSE IF ((r s (z w ) > p s AND r s (zi ow ) -p s <a c AND a c < r s (zi ow ) + p s ) 
OR (p s > r s (z hw ) AND p s - r s (zi ow ) < a c AND a c < p s + r s (z, ow ))) 

THEN ! Subcase 4.4 

z high = MIN[z, + yjaj - (a c - p s f , z,] 

IF (a s >a c + p s ) THEN 

z hi g h 



V in =V in + jdzA in =V in + AV 3 VIIIA.4(9) £176) 



where AV3 is given in the next section below. 
ELSE 



thigh 

V in =V in + \dzA in =V in + AV 4 Vm.A.1(l0) O77) 



where AV4 is given in the next section below. 
END IF 

IF {zhigh = z t ) keepGoing = FALSE 

END IF ! End Subcase 4.4 

END IF ! End Subcases 

END IF ! End hemisphere IF's 

ziow - Zhigh ! Prepare for next z interval 

END DO ! End DO WHILE loop 

Z high 

Evaluation of the Integral Evaluation of the Integral J dzA in 

Z low 

[00378] We now address the integral 
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J dzAj, 



)078} 



where A in is given b yEqo. VIII.A. 4 (2d, 3, and 4 ) Eqs. 169d. 170a, 170b. 171a. and 171b as 



A in = al cos" 1 



2psa c 



+ (af - £ 2 ) cos 



-i 



f p] - fl » + fl»-^ 



\ AC + (A + O 2 - «. a ][« s 2 " (A " a,) 2 " C] 



with 



)(179} 



)08Q) 



[00379] This integral will involve elliptic integrals of the first, second, and third kinds 
which for present purposes are defined as in Abramowitz and Stegun. The incomplete elliptic 
integral of the first kind is given by 



sin 2 0 



)am 



[00380] The complete elliptic integral of the first kind is denoted by K(k) and is equal to 
F(nl2,k). Similarly, the incomplete elliptic integral of the second kind is given by 



E(<p,k) = \dO y]l-k 2 sm 2 0 



with the complete elliptic integral of the second kind denoted by E(k) and equal to E{nl2,k). 
Finally, the incomplete elliptic integral of the third kind is given by 
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Tl(<p,n,k)= $d0 



(l-nsm 2 0)yjl-k 2 sin 2 0 



VIII.A.K16 X 183) 



with the complete elliptic integral of the third kind denoted by H(n,k) and equal to Yl(ii/2,n,k). 
Various efficient routines for the evaluation of these special functions can be found in the 
literature. 

[00381] Focusing on the integral of each of the three terms in Eq. VIIL A. 4 (12), 179, it 
is possible to perform an integration by parts of the first term of Eq. VIII.A. 4 (12) 179 to 
dispatch the inverse cosine term. 



<T2 

a; jac cos 

^2 £ 2 



2p s a c 



f ^ 



= a c 2 ^cos 



-1 



ps + a 2 c -a* +C 
2p s a c 



2^ 



)084) 



where 



R U2 = y/[a? - (a c - p,f - C 2 ][C+ (a e + p s f - a 2 ] 



The right-most integral of Eq. VIH.A. 4 (17) 184 is an elliptic integral. 

[00382] An integration by parts also may be used for beginning the integration of the 
second term of Eq. vm.A.4fl2V 179. 
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f2 



f p) -a] + a 2 -C^ 



\ 2p s 4a}-C 



J 



f 



2\ 



' 3 



COS" 



r p]_-al+a^^ 
2p s yja? -C 



2\ nl/2 



3(a, 2 -0*' 



<T 2 



a. 2 - — cos 1 



2 -a 2 +a f 2 
2p^a]-C 



VIII.A. 4 (19) (186) 



- (a 2 + 3a, 2 - p) )\d£ 4tt + - %C " - ^ 3 (flc - A 2 / f d£ — — - 
3 V H "l b R}' 2 3 1 h R U2 3 (fit -OR 



2\ Dl/2 



where the integrals after the second equal sign in Eq. VIII.A. 4 (19') 186 come from partial 
fraction separation of the integrand after the first equal sign. Each of the integrals involves one 
or more elliptic integrals. 

[00383] The third term of Eq. VIII.A. 4 (12) 179 is already in a form amenable to 
solution by elliptic integrals: 



4 K* 1/2 • 



)£187) 



[00384] Combining these three expressions in Eqs. VIII.A.<1fl7. 19. and 20). 184. 186. 
and 187 y ields the general expression for the A Fused in Eqs. VIII.A. 4 (7 to 10): 174 to 177: 



AF = a 2 ^ 2 )cos- 



f ~2 



pj + a 2 - a? +£ 



2\ 



cos 



2p 5 a c 

\ f 2 



(2 



f 1 ^ 

2 1 n 



V 3 



J 



+ -a?(a? -p 2 )L 0 + 



|a 2 -a} + I A 0z 2 + Il 4 - X - a j{a 2 c - p j)P-i E 



mm 



where 
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021} 




and 

E = \d<ZR xn Vm.A. 4 f26) fl93) 

with R m given by Eq. VIII.A. 4 (18V 185. The explicit solutions to the L 9 P, and E integrals of 
Eqs. Vm.A.4(22 26) 189 to 193 depend on the values of C/, G, a S9 and (a c - p s ). These 
expressions are used in the preferred implementation. 

Explicit Solution for Interior Volume When a 5 ->-g e — -p JExplicit solution for 
interior volume Vm when a z > a c - p 1 

[00385] In the case where a s 2 > (a c + p s ) 2 , the preferred implementation follows the (* 
integration up the z axis, the southern hemisphere's sphere circle eventually circumscribes the 
cylinder circle unless it hits the top of the cylinder first. The limits of integration are therefore 
given by 

C, = MAX[B 9 -(z t -z s )] and ( 2 = MJN[A, -(z b - z s )] VHLA.4(27) (194) 
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where A and B are given below. If a s 2 > (a c + p s ) 2 , then it follows that a s 2 > (a c - p s ) 2 . Therefore, 
we define 



A = yja] -(a c -p 5 ) 2 



)£195a) 



and 



B = yja] -(a c + p s ) 2 



)£195b) 



and note that A>B>0. Then 



R" 2 = J(Ai-C>)(C-B>) 



[00386] The integrals for this form of R m are readily found in references such as 
Gradshteyn and Ryzhik. Define 



\A 2 -B 2 



)097a} 



h = sin I ^- L - 

VA 2 -B 2 



)(197b) 



and 



Ja^b 2 

q= —r- 



)(197c) 
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[00387] Then, recalling the definitions for the elliptic integrals of Eqs. VIII.A. 4 (1 4 to 
i€h - 181 to 183, t he following expressions for the integrals are obtained: 



L 0 = A A [Y{X x ,q)-E{X 2 ,q)l £><A 



)098) 



Note that E(X 2 ,q) = 0 if £> = A. 



L 2 =A[E{X x ,q)-E(X 2 ,q)l 6<A 



)£199) 



Again, E(X 2 ,q) = 0 if £> = A. 



U = A/3 [2(A 2 + B 2 ) E(X,,q) - B 2 ¥(X h q)] + fy/3 [(A 2 - Cj 2 )(C, 2 - B 2 )] m - 
A/3 [2(A 2 + B 2 ) E(X 2 ,q) - B 2 F(X 2 ,q)] + 6/3 [(A 2 - C 2 2 )(C 2 2 - B 2 )] 112 , £ 2 <A VIII.A.4(33) (200) 



Note that the entire second line is zero if C2 = A. 



P = 



1 



A{a 2 -A 2 ) 



n 



- A 2 -B 2 
K— r>9 



^ ( 



A 2 -- 2 



a 



n 



A 2 -B 2 



A 2 , 



A 2 -a 2 



i j 



)O0_u 



The second elliptic integral-wittbk - with X? as its first argument is zero \iQi=A. 
[00388] Lastly, 

E = A/3 [(A 2 + B 2 ) E(X h q) - IB 2 E{X,,q)] - C//3 [{A 2 - tfxtf - B 2 )] 1 ' 2 - 

A/3 [(A 2 + B 2 ) E(X 2 ,q) - 2B 2 F(X 2 ,q)] - 6/3 [{A 2 - C 2 2 )(6 2 - B 2 )} 112 , Q 2 < ,4 VIIIA. 1 f35¥ 202) 
Again, the entire second line of this expression is zero if ( 2 = A. 
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Explicit Solution for Interior Volume V j » When a s^rUe^-frs Explicit solution for 
interior volume Vi n when a 1 <a (1 + p 1 

[00389] In the case where a 2 < (a c + p s ) 2 , A 2 is still defined as it was before, A 2 = a 2 - 
{a c - Psf> but B 2 becomes the negative of its previous definition: 



B 2 = (a c +p s ) 2 -a s 2 . 



¥203) 



[00390] Although A 2 is always greater than zero, we simply have that B 2 > 0. 
The radical is defined differently as 



R U2 = y j(A 2 -£ 2 )(C 2 + B 2 ) 



(204) 



and that changes the elliptic integral solutions of Lq, L2, L 4 , P, and E. Defining 



y x = sin " 



A \B 2 + £\ } 



){205a} 



Yi = sin 



Ci A 2 + B 2 



A VB 2 + C 2 



2 J 



){205b} 



r = 



4a t +b 2 



){205c) 



the following expressions are obtained for the integrals. Note that & is inside absolute value 
signs so that only its magnitude is used. The upper limit however, may be positive, negative, 
or zero. A negative & results in a change in sign for all terms in the integrals involving All 
the integrands are even in C so that all the integrals are odd. The elliptic integrals are also odd in 
their first argument. The expressions for the desired integrals are now 
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L 0 = 



[F(y x ,r) + F{y 2 ,r)} 



){206J 



L 2 = jA 2 + B 2 E{y x ,r)- 



B 2 „, . \A 2 -t;\ 



F(y u r)-\^\ 



B 2 + a 



JA 2 + B 2 E(y 2 ,r)- 



B 2 „. . „ M 2 -£ 2 



){207} 



I 4 = 



3^ 2 + 5 2 



[(22? 2 - A 2 )B 2 F(y v r) - 2(B* - A*)E{y„r)}- 



a' -a 

2 . /-2 



+ 



1 



Vm.A.4(4 1 X 208) 



3y[A 2 + B 2 



[(2B 2 - A 2 )B 2 F(y 2 , r) - 2(B* - A<)E(y 2 , r)] - 



&(2A 2 -B 2 + £ 2 )l^—£- 
3 2 Vfl 2 + £ 



P = 



a 2 (a] +B 2 )JA 2 + B- 



:[B 2 Yl(y x ,n,r) + a) F(y x ,r) + B 2 Yl(y 2 ,n,r) + a 2 F(y 2 ,r)Y- 



)Q09a) 



where 



n = 



A\a) +B 2 ) 
a 2 (A 2 + B 2 ) 



)(209b} 



and 



E = ^A 2 + B 2 [B 2 F(y„ r) - (B 2 - A 2 )E{ 7l , r)] + ^ (tf + IB 2 - A 2 \f^^ 
3 3 \B 2 +Q 

+ U A 2 + B 2 [B 2 F(y 2 ,r)-(B 2 - A 2 )E(y 2 ,r)] + + 2B 2 - A\l££- 

3 3 y B 2 + £ 2 
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which completes the solutions for the intersection of a sphere and cylinder. The explicit 
expressions for & and & ar e given in S e ction VIIIAAa. earlier. These expressions are used in 
the preferred implementation to determine the intersection. 

[00391] The task of determining the packing fraction also may be accomplished 
according to an aspect of the invention using a modified approach. 

[00392] Recall that the radial distance of the sphere center from the z axis by 



[00393] Also recall that there are four possible cases for this overlap problem. In 
Case 1, the sphere is completely outside the cylinder so V in = 0. In Case 2, the sphere is 
completely inside the cylinder so V in = (4/3 )na s 3 . In Case 3, the cylinder is completely inside the 
sphere so V in = nac 2 (z t - z b ). In Case 4, the sphere and cylinder intersect each other and the 
expression for V in is more complex as will be derived below. 

Case 1: — Sphere Completely Outside Cv/*;i<fc r Case 1: Sphere Completely Outside 
Cylinder 

[00394] The sphere is completely outside the cylinder 



Ps = [x s 2 +y s 2 ] 



m-Z) (2U) 



0 



03-5 ¥212) 



only if one of the following conditions are met: 



If the sphere's north pole is lower than the cylinder bottom . . . 



z 5 + a s < z b 




Or if the sphere's south pole is above the cylinder top . . . 
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z s -a s >z t (2.3 6b Y 213b) 

Or if the sphere is more than its radius away from the cylinder side . . . 

p s -a s > a c (2.3 6c) (213c) 

Or if the sphere center is above the cylinder top and is more than its radius away 
from the upper cylinder edge . . . 

ips - o. c f + (z s - z t f > a s 2 and z s >z t and p s >a c (2.3 6d) (213d) 

Or if the sphere center is below the cylinder bottom and is more than its radius 
away from the lower cylinder edge . . . 

„ 2 



(ps-a-c) +(z b -z s ) >a s and z s <z b and p s >a c . (2.3 6 e ) (213e) 

Case 2: Sphere Completely Inside Cv/i/irfgr Case 2: Sphere Completely Inside Cylinder 
[00395] The sphere is completely inside the cylinder 

^•„ = (4/3)W 0^-3¥214) 

only if all three of the following conditions are simultaneously met: 

If the sphere's north pole is no higher than the cylinder top . . . 

z s + a s < z t (2.3 4 a) (215a) 

And if the sphere's south pole is no lower than the cylinder bottom . . . 

z s -a s >z b (2.3 4b) (215b) 
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And if the sphere's equator nowhere exceeds the cylinder radius . 



ps + a s <a c • (2.3 4 c) (215c) 

Case 3: Cylinder Completely Inside Spher e Case 3: Cylinder Completely Inside Sphere 
[00396] The cylinder is completely inside the sphere 



Vin = na 2 (z t - z b ) (2^-7X216) 



if both the following two conditions are met: 

If the upper edge of the cylinder is inside the sphere . . . 

(p s + a c f + (z, - z s f < a 2 (2.3 8a) (217a) 

And if the lower edge of the cylinder is inside the sphere . . . 

(p s + a c ) 2 + (z s -z b ) 2 <a s 2 . (2.3 8b) (217b) 

Case 4: Cylinder and sphere interscc t Czse 4: Cylinder and Sphere Intersect 

[00397] Again, there are several subsets of Case 4. To efficiently organize them, we 

define two circles: the sphere circle and the cylinder circle. The sphere circle is the perimeter of 

the sphere's cross section at a given z altitude. Its radius is 



r s = [a 2 -(z-z s ) 2 ] U2 , z s -a s <z<z s + a s ft^-9 ¥218) 



and its radial position p s is given by Eq.42V 211. Likewise the cylinder circle is the perimeter of 
the cylinder's cross section which has constant radius a c at any z altitude between z b and z t and is 
centered on the z axis. 
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[00398] To describe the intersection of the sphere and cylinder, we define four subcases 
fe^he- for the relationship of the sphere circle and the cylinder circle at a particular altitude z: 
Subcase 4. 1 : The sphere circle is completely enclosed in the cylinder circle. 
Subcase 4.2: The two circles do not overlap or one or both are nonexistent. 
Subcase 4.3: The cylinder circle is completely enclosed in the sphere circle. 
Subcase 4.4: The two circles intersect. 
[00399] To find our task is reduced to integrating in z the overlap areas of the two 
circles A in from the altitude where one subcase holds to the altitude where another subcase holds. 

For Subcase 4. 1 , A in = na s 2 . (2.3 - 10a) (219a) 



For Subcase 4.2, A in = 0. (2.3 10b) (219b) 

For Subcase 4.3, A in = na c 2 . (23 10c) (219c) 

For Subcase 4.4, A in = a c (a- cosa sina) + r s (fi - cos/? sin/?) 




where r s is given by Eq. (2 3 9) 218 and 



2 2 2 

cos« = A +<Xc ~ r * (2.3 Hal f 220a) 



q 2 j_ v 2 — or 1 

cos p = ti 1 s- (2.3 12a) (221a) 



2A^ 



[00400] It is now necessary to determine a series of altitudes at which the cylinder circle 
and the sphere circle change from one subcase to another in Case 4. 
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[00401] A sphere and a cylinder side can intersect each other at either 0, 2, or 4 altitudes. 
Pairs of altitudes may be degenerate (two solutions having the same altitude) if the sphere just 
nicks the cylinder wall. The solutions for these altitudes come from solving the following 
equation: 

[r s (z)f = [oc ± p s f (2.3 13a) (222a) 

or 

a s 2 - (z - z s f = [a c ± p s f (2.3 13b) (222b) 

which, using Eq. (2 3 9), 218, has the four solutions 

z = z s ±ja> -(a c ±p s ) 2 . (2^-44)Q23} 

The + sign in the first ± of Eq. (2 3 1 4 ) 223 denotes the upper or northern hemisphere 
intersection altitude; the - sign denotes the lower or southern hemisphere intersection altitude. 
The + sign in the second ± of Eq. (2 3 1 4 ) 223 denotes intersection with the far cylinder wall, 
on the opposite side of the z axis from the sphere's side; the - sign denotes intersection with the 
near cylinder wall, on the same side of the z axis as the sphere. 

[00402] Armed with these four solutions, we will now consider each of the four 
subcases of Case 4, but we will consider separately whether the sphere circle lies in the southern 
hemisphere of the sphere or whether it is the equatorial plane or northern hemisphere. We 
consider the four subcases separately for the two hemispheres making eight situations in all. We 
consider the hemispheres separately because if the sphere circle is in the southern hemisphere, it 
is increasing in size as z increases. Hence the preferred implementation can determine the 
altitude at which subcases change. If the sphere circle is its equatorial plane o r in its th e is its 
northern hemisphere, it is decreasing in size as z decreases, and the altitude at which subcases 
will change again can be determined. 
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[00403] We now outline the overall logic used in the preferred implementation to find 
the altitudes at which the subcases for Case 4 change. Let z\ ow be the lower altitude of the 
subcase region currently being examined and Zhi g h be its upper altitude. 

[00404] The absolute bounds boundary below o r abov e which th e r e can above, with 
neve r b e any any overlap between sphere and cylind e r cylinder, is given by 

z 0 = MAX[z b , z s - a s ] (23 151 (224) 

for the lowe r bound boundary and 

z 5 = MIN[z,, z S9 + a s ] (23 16 X 225) 

for the uppe r bound, boundary. The altitudes z\ through Z4 correspond to the different 
intersection altitudes (if they exist) of the sphere's polar plane that is perpendicular to the 
cylinder axis. A sphere's polar plane is the plane which intersects the sphere in a great circle 
perpendicular to the equatorial plane. It can be oriented at any azimuthal angle. The polar plane 
of interest to us is the one passing through the cylinder axis. Th e int e rs e ction altitud e s ar e 
illustrat e d in Fig. - 

[00405] In ascending order (again, if they exist) the intersection altitudes are: 

Southern hemisphere intersection with cylinder wall closest to sphere-center 
center: 

z x = z,-ylat -(a c ~p s f (2^4^(226} 

Southern hemisphere intersection with cylinder wall furthest from sphere-eenter 
center: 

*2 = *,-A 2 "( fl c + A) 2 (2.3 1S)( 227) 
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Northern hemisphere intersection with cylinder wall furthest from sphere-center 
center: 



Northern hemisphere intersection with cylinder wall closest to sphere c e nt e r 
center: 

z 4 = z s + y Jal -{a c -p s f . 

In all four of these solutions, if the argument of the square root is negative, the square root is 
replaced by zero. If the sphere intersects both sides of the cylinder, there are four real solutions; 
if only one side, there are two real solutions; if not at all, there are no real solutions of Eqs.-I^? 
through 20. 226 through 229. The figure below shows the two real solutions for a sphere that 
intersects only one side of the cylinder. 

[00406] The method used in the preferred implementation for finding the volume V in of 
a sphere that lies inside a cylinder will now be discussed. The method is divided into two DO 
WHILE loops, one for the southern hemisphere and one for the northern hemisphere (which 
includes the sphere's equatorial plane). This is convenient because for the southern hemisphere 
the sphere circle is always increasing with increasing z while for the northern hemisphere it is 
always decreasing. The method as specifically implemented here has a number of comments so 
that the reader can follow the reasoning guiding it. 

Vi n = 0 ! Vj n is the volume of this sphere lying inside the 

! cylinder. Initialize it to zero, 
zo = MAX[zb, z s - as] ! zo and z$ are the absolute lower and upper 

z 5 = MIN[z t , z s + as] ! limits of possible sphere-cylinder overlap. 

EF (z s > Zb) THEN ! The southern hemisphere can contribute to Vj n if and 

! only if this statement is true. Otherwise it lies 
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! beneath the cylinder bottom and cannot contribute. 
! (It might not contribute even when the above 
! statement is true but that is checked below.) 

ziow = zo ! Initialize zi ow to smallest possible value. 

r s (ziow) = [as 2 -(ziow-z s ) 2 ] 1/2 ! Find initial sphere circle 

! radius. It is 0 if square 
! root argument < 0. 

! Find the initial southern hemisphere subcase. 
IF (p s + r s (ziow) < ac) THEN 

subcase = 1 ! Sphere circle inside cylinder circle. 
ELSE IF (p s - r s (z low ) > ac) THEN 

subcase = 2 ! Sphere circle outside cylinder circle. 
ELSE IF (r s (zi ow ) - p s > ac) THEN 

subcase = 3 ! Cylinder circle inside sphere circle. 

ELSE 

subcase = 4 ! Sphere circle intersects cylinder circle. 
END IF 

! We exit the DO WHILE loop when the subcase is no longer a positive 
! integer. Remember that since we are in the southern hemisphere, the 
! sphere circles are expanding with increasing z. 

DO WHILE (subcase > 0) ! Southern hemisphere loop. 

IF (subcase = 1) THEN 

! The sphere circle lies inside the cylinder circle. 
! It expands with increasing z until one of three 
! things happen: 

! (1) It hits the top of the cylinder. 
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! (2) It intersects the cylinder circle. 

! (3) It reaches its equatorial plane. 

! Numbers (2) and (3) can be combined if we agree to 

! define z\ with the square root replaced by zero if 

! its argument is negative. 

w = as 2 - (ac - p s ) 2 
IF (w > 0) THEN 
zi =z s - w 1/2 

ELSE 

Zi = z s 
END IF 

Zhigh = MIN[z t , zi] 

+ Vi(zi 0W ,Zhigh) ! We'll derive an explicit 
! expression for this Vj later. 



! Now we must find the next subcase. 

IF (z high < MIN[z s ,Zt]) THEN ! There is another 

! subcase to be 
! addressed in the 
! southern hemisphere. 
IF (p s > 0) THEN ! Sphere is not concentric with 
subcase = 4 ! the cylinder so it will 

! expand until it intersects. 
ELSE ! Here, the sphere is concen- 

subcase - 3 ! trie so it will expand 
END IF ! until it perfectly circum- 

! scribes the cylinder 
! circle. 

ELSE ! You reached the end of southern 
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subcase = 0 ! hemisphere contribution. End the 
END IF ! S. hemisphere's DO WHILE loop. 



ELSE IF (subcase = 2) THEN 

! The sphere circle lies outside the cylinder circle 

! until it expands to the point where it intersects 

! the cylinder. This intersection must happen or we 

! wouldn't be here in Case 4. The intersection 

! happens at z\. There is no contribution to Vi n from 

! this region. We only need to find Zhigh = Zj. 

w = as 2 - (ac - p s ) 2 
IF (w > 0) THEN 

Zhigh - z s - w 

subcase = 4 ! Not done yet with S. hemisphere. 

ELSE 

Zhigh ~ Z s 

subcase = 0 ! Reached equatorial plane. Now 
END IF ! we're done with S. hemisphere. 

ELSE IF (subcase = 3) THEN 

! The cylinder circle lies inside the sphere circle. 
! Since the sphere circle is expanding, it will 
! remain inside until it reaches the first of either 
! the top of the cylinder or its equatorial plane. 

Zhigh = MIN[z t , z s ] 

Vin = Vi n + V2(zi 0W >Zhigh) ! We'll derive an explicit 

! expression for this V2 later. 

subcase = 0 
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ELSE IF (subcase = 4) THEN 

! The only possibility left is that sphere and 

! cylinder circles initially intersect. This region 

! extends upward until one of three things happen: 

! (1) We reach the top of the cylinder. 

! (2) We reach the sphere's equatorial plane. 

! (3) We reach a point at which the sphere circle 

! completely circumscribes the cylinder circle. 

! We cannot reach the top of the sphere in this case 

! because we are only dealing with the southern 

! hemisphere. (Remember z 2 = z s if the argument of 

! the square root is negative.) 

w - a s 2 - (ac + p s ) 2 
IF (w > 0) THEN 

z 2 = z s - w 1/2 

subcase = 3 

ELSE 

Z2 = Z S 

subcase = 0 
END IF 

Zhigh = MIN[z t , z s , z 2 ] 

Vj n = Vj n + V3(zi ow >Zhigh) ! We'll derive a numerical 

! integration for this V3 later. 

END IF 

! If you have reached either the cylinder top or the 

! sphere's equatorial plane, you are done with the southern 

! hemisphere's contribution to V in . Go to the northern 
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hemisphere next. (We use > signs below in case of round- 
off error. With infinite precision, Zhi g h could never be 
greater than z t or z s .) 



IF (Zhigh > z t OR Zhigh > z s ) THEN ! Redundant check saying 
subcase = 0 ! we're done with the southern hemisphere. 

ELSE 

ziow = Zhigh ! Cycle back through the DO WHILE loop 
END IF ! for the next region in the southern 

! hemisphere contributing to Vj„. 
END DO ! End of southern hemisphere DO WHILE loop. 



END IF ! End of southern hemisphere overall EF statement. 



IF (z s < z t ) THEN ! The northern hemisphere can contribute to Vi n if and 

! only if this statement is true. Otherwise it lies 

! above the cylinder top and cannot contribute. (It 

! might not contribute even when the above statement 

! is true.) 



ziow = MAX[z b , z s ] 

rs(ziow)=[as 2 -(ziow-Zs) 2 ] 1/2 ! Find initial sphere circle radius in 

! N. hemisphere. It is 0 if square root 

! argument < 0. 



! Find the initial northern hemisphere subcase. 
IF (p s + r s (z low ) < ac) THEN 

subcase = 1 ! Sphere circle inside cylinder circle. 
ELSE IF (p s - r s (z, ow ) > THEN 

subcase = 0 ! Sphere circle outside cylinder circle. 
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! Although this is actually subcase 2, since 

! N. hemisphere sphere circle can only contract 

! it will never intersect the cylinder 

! circle if it does not initially intersect it. 

! Therefore, we set the subcase to zero so the 

! N. hemisphere's DO WHILE loop will end. 
ELSE IF (p s + ac < r s (z ]ow )) THEN 

subcase = 3 ! Cylinder circle inside sphere circle. 

ELSE 

subcase = 4 ! Sphere circle intersects cylinder circle. 
END IF 



DO WHILE (subcase > 0) ! Begin northern hemisphere loop. 



IF (subcase =1) THEN 

! The sphere circle lies inside the cylinder circle. 
! Since the sphere circle shrinks as z increases, it 
! cannot ever intersect the cylinder circle. It can 
! only end at its north pole or at the top of the 
! cylinder. That ending position is Z5. 



Zhigh - Z5 

+ Vi(ziow>Zhigh) ! We'll derive an explicit 
! expression for Vi later. 

subcase = 0 



ELSE IF (subcase = 2) THEN 

! The sphere circle is inside the cylinder circle. 
! In the northern hemisphere this means they will 
! never intersect more so there will be no more 
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! contribution to Vj n . 
subcase = 0 

ELSE IF (subcase = 3) THEN 

! The cylinder circle lies inside the sphere circle. 
! The sphere circle will contract until we hit top 
! of the cylinder or until we reach Z3 which is the 

! point at which the two circles will intersect. 



w = a s 2 - (ac + p s ) 2 



! w is nonnegativ e non-negative or we 
! couldn't be here. 

z 3 = z s + w 1/2 
Zhigh = MIN[z t , z 3 ] 

+ V 2 (zio W? Zhigh) ! We'll derive an explicit 
! expression for V 2 later. 

IF (z hign <z t ) THEN 



IF (p s > 0) THEN 
subcase = 4 



ELSE 



subcase = 1 



END IF 



ELSE 



subcase = 0 



END IF 



! Sphere is not concentric 
! with cylinder so it will 
! contract until it inter- 
! sects the cylinder circle. 
! The sphere is concentric so 
! it will contract until it 
! suddenly lies completely 

! inside the cylinder circle. 
! If we're here, we reached 
! top of the cylinder; we're 
! done with N. hemisphere. 



ELSE IF (subcase - 4) THEN 
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! The only possibility left is that sphere and 
! cylinder circles initially intersect. This region 
! extends upward until one of two things happen: 
! (1) We reach the top of the cylinder before 
! reaching the north pole. 
! (2) The intersection stops with the sphere circle 
! now either completely inside or completely 
! outside the cylinder circle. This happens at 
! z 4 . 

w = as - (ac - p s ) ! w is never negative or we 

! wouldn't be here. 

z 4 = z s + w 1/2 
z high = MIN[z t , z 4 ] 

Vj n = Vin + V 3 (ziow ? Zhi g h) ! We'll derive a numerical 

! solution for V 3 later. 

IF (p s < ac) THEN 

subcase = 1 ! Next, sphere circle will lie 
! inside cylinder circle. 

ELSE 

subcase = 0 ! Next, sphere circle will lie 
! outside cylinder circle. 

END EF 
END IF 

! (Note the second logical expression in the IF statement below is z 5 , 
! notz s .) 

IF (Zhigh > z t OR Zhigh > z 5 ) THEN ! Redundant check saying 
subcase = 0 ! we're done with northern hemisphere. 

ELSE 
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ziow = zwgh ! Cycle back through the DO WHILE loop 
END IF ! for the next region in the northern 
! hemisphere's contribution to Vi„. 

END DO ! End of northern hemisphere DO WHILE loop. 
END IF ! End of northern hemisphere overall IF statement. 

[00407] Expressions for the three volume expressions F,(z/ 0H , 5 z/,^), i = 1, 2, or 3 will 
now be addressed. Remember in the following that a c is the cylinder radius, a s is the sphere 

radius, (x S9 y S9 z s ) is the position of the sphere's center, and p s = <Jx 2 + y 2 s is the radial distance 

of the sphere center from the cylinder axis. 

Derivation of Vg JDERIVATION OF V L 

[00408] The volume V\ is an integral of full sphere circles from z/ ow to z high because the 
sphere circles lie entirely within the cylinder circle. 

z high z high 

V, = n \dz[r s (z)f = x \dz[a 2 5 -(z-z s ) 2 ] Q^gm 

z low z fow 

or 

^i=y(^-^J[3^ 2 -z 2 high -zf ow + z high (3z s - z lo J - 3z s (z s - z lo J] (2.3 22)( 231) 

where z/ ovv and Zhigh depend on where we are in the above method. Their values are given in the 
implementation shown above. 

Derivation of V ^DERIVATION OF V z 

[00409] The volume V 2 is the integral of the cylinder circle from z iow to z h i g h because the 
cylinder circle lies entirely within the sphere circle. It is given simply as 
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V 2 = 7ia\ (z high - z low ) 
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Derivation of V ^ DERIVATION OF V, 

[00410] The volume Vy is the integral of intersecting sphere and cylinder circles from 
zi ow to Zhi g h. The limits are always given in the method as described above. The intersect area is 
given by Eqs. 2 3 10, 11. and 12 219, 220, and 221 which we must integrate in z. 



'high 

F 3 = J dz{a] [a - cos a sin a] + [af - (z - z s ) 2 }[/3 - cos fi sin /?]} (2.3 2 4 V 233) 



where 



cos a = 



2PA 



)£234) 



sin a = 



_ > 2 -te-A) 2 -(*-*,) 2 ][fa + A) 2 -fl, 2 +(g-g.) 5 ] 



){235} 



a = arccos 



Ps + <*! 



2Ps"c 



(2.3 27X 236) 



and 



cos P = 



Ps -a 2 +Q 2 -(^-^) 2 
2p> 2 -(z-z,) 2 



)Q37) 



sin /? = 



Vfo 2 - (a c - p 5 f - (z - zj 2 ][(a c + A ) 2 + a] - (z - zj] 
2aV« 2 -(*-*.) 2 



)Q38) 
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P = arccos 



2p^a) -{z-z s f 



(2.3 30 ¥ 239) 



[00411] The integral of Eq. 2 3 2 4 233 could be done analytically in terms of complete 
and incomplete elliptic integrals of the first, second, and third kinds. However, this solution is so 
cumbersome that it is more efficient to perform it numerically, especially because the integrand 
is quite smooth and well-behaved so that a Romberg integration is very efficient and contains an 
accurate error assessment. 



Calculating the Total Volume Fraction for a Paeft Calculating the Total Volume 
Fraction for a Pack 

[00412] The preferred implementation also can determine the volume fraction c|> of the 
pack. Only the innermost cylinder contains particles that look like a full three-dimensional pack. 
The outer cylinders do not contain the particles the inner cylinders contain. However, it would 
be unfortunate to throw away all the pack information the outer cylinders contain. They do 
perturb the pack of the innermost cylinder correctly and their correct geometry for the particles 
they contain should somehow be used in estimating the volume fraction of the pack. 

[00413] The preferred implementation determines the volume fraction ^ in the 
following way. The y-th cylinder can contain all particle modes whose cylinder radii are less 
than or equal to its cylinder radius a c (j). As long as the water level option of pack building is 
operating, the y-th mode particles are distributed approximately as they would be in a truly 
three-dimensional pack. The preferred implementation therefore calculates the volume fraction 
<f> } of Mode j particles within the y-th cylinder and calls it the y-th partial volume fraction. It does 

likewise for all the other modes, e.g., J of them. Then the total volume fraction is simply the 
sum of all the partial volume fractions: 




If desired, the user can specify that only the volume fraction of the innermost cylinder, where the 
pack is three dimensional, will be calculated. 

[00414] To calculate the partial volume fraction <j> } in the y-th cylinder, the preferred 

implementation sums the interior volume of each Mode j sphere inside the y-th cylinder. Let nj 
be the number of Mode j spheres in the y'-th cylinder. Let Vj(k) be the volume inside the cylinder 
of the &-th sphere of Modey. Then 

Ui) = iw • VIII.A.5(2) (241) 

k=\ 

[00415] Dividing this total interior partial volume by the volume of the y-th cylinder 
V cy i{j) gives the partial volume fraction <j>j for the pack: 

<f>j =V tot {j)IV cyl (j) VIIIA.5(3) (242) 

where 

V C yi(j) = n[a c (j)] 2 [z top (j) - z bot {j)\ VIIIA.5f4) (243) 

[00416] The user has the option of making calculating volume fractions in cylinders that 
are a little smaller than the actual cylinders used in making the pack so that possible wall effects 
can be eliminated from the calculation. Remember, however, that the cylinder walls are barriers 
to the particle centers, not their leading edges so that disturbances to the pack from wall effects 
normally will be small. 

Finding the Radial Distribution Functions Findinz the Radial Distribution Functions 
[00417] The radial distribution functions, denoted by gy(r)dr, are probability distribution 
functions denoting the mean number of Mode or Submode j particles whose centers are between 
r and r + dr from the centers of Mode or Submode / particles throughout the pack. They can be 
important in establishing many properties of random packs. Because the packs here are assumed 
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random, only the scalar radial distance r is used rather than vector r in this implementation. 

[00418] Reduced dimension packs do not extend equally in all directions. They are a 
series of concentric cylinders. Therefore, for a given Mode or Submode / particle, the preferred 
implementation calculates the volume of the spherical shell segment centered on Particle i and 
lying between r and r + dr, that lies inside the y-th cylinder. Then it counts all Mode j particles 
lying in that shell segment. Dividing the number of Mode j particles by the volume of the shell 
segment gives an approximation to gij(r)dr at that position r. 

[00419] Finding the volume of this spherical shell segment is not a trivial procedure. 
Fortunately, a means to calculate the total sphere volume lying inside any cylinder is described 
above. Therefore, the preferred implementation can calculate the total volume of an imaginary 
sphere of radius r + dr, centered at Particle /, that lies inside the 7-th cylinder, and subtract from it 
the total volume of the concentric sphere of radius r to obtain the volume of the spherical shell 
segment between r and r + dr r + dr lying in the y'-th cylinder. 

[00420] The user specifies the maximum radial distance r max {ij) for which gy is desired 
and the number of bins nty into which it divided r max (ij). Then the bin widths have finite width 

Ar /y = r max {ij)lmij Vm.B(l) (244) 

rather than infinitesimal width dr. Then the g tj are determined as follows. 

[00421] First, note that the gij{r) are zero for r < a t ■ + cc,. Then, consider a particular 
Mode / particle with index number i h . Denote the number of Mode j particles in the first bin 
from Particle k particle by NJJ, h \ a ( - + aj). The first bin lies between (cc, + aj) and (a, + aj) + Ar,y. 
In general, the number of Mode j particles lying in the &-th bin from Particle ih is denoted by 
Nj(i h \ r k ) where r k = a ( + aj + (£-l)Ar £ y, k = 1, 2, . . . , my. 

[00422] Now consider the volume of the k-th bin from Particle ih that lies inside the j-th 
cylinder. This can be obtained from the volumes V in of previous sections where V in = V in (x S9 y s , 
z s > a s \ «c, Zbot, Ztop)- The argument list contains the position of the sphere center (x Si y s , z s ) 9 the 
sphere radius a Si the cylinder radius a c? and the cylinder's bottom and top altitudes (zb oh z top ). 
The volume for the k-th bin about Particle ih is then given by 
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Vintthj, k) = V in (x ih , y jh , z ih , r k +r, a cj , z botJ , z lopJ ) - V in (x ih ,y ih , z ih , n; a cJ , z bolJ , z MrJ ) Wlll.B(2) (245) 



and the number density of Mode j particles in this volume is 

<r(ihj; k) = Nj{i h ; r k )IV in {i h J; k) . VIII.B(3 ¥ 246) 

[00423] Finally, to obtain the gi/r), the above expression can be averaged over all Mode 
/ particles: 



gij (r) = - Jo-fey; k) VIII.B( 4 )( 247} 



if r > a, + Oj and is zero otherwise. The term k is the integer lying in the range 

(r - a t • aj)/Anj < k < 1 + (r - a, - tf/)/Ar /y . VIILB(5) (248) 

If r is so large or so small for a particular particle //, that the spherical shell containing it misses 
the y'-th cylinder altogether, then that term in the sum of Eq. VIII.B( 4 ) 247 is omitted altogether. 
[00424] The normalized radial distribution function, denoted G^r), is given by 



Gi/r) = gij(r)/(47tr 2 ) VIII.B(6) (249) 



or in our finite bin width treatment, it is given by 



Gfr) = 4 & " K J . VIII.B(7)Q 50} 

-<r 3 k + x -rl) 



It is the mean number density n } of Type j particles at a distance r from a Type i particle. Thus 
normalized, the function G,y(r) approaches the constant «, as r grows large rather than increase 
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quadratically as gy(r) does. -*y-j^_represents the number-off - of / p articles per unit volume, e.g., 
concentration, output from the preferred implementation. 

Finding the Coordination Numbcrs Y'mAmz the Coordination Numbers 
[00425] The coordination number Qy is the mean number of Mode j particles touching a 
Mode / particle. Again, because a reduced dimension pack is ordered in concentric cylinders 
rather than being isotropic in all directions, the preferred implementation finds how much of the 
surface of the z-th particle is available for Mode j particles to touch. That is equivalent to asking 
how much of the imaginary spherical surface centered on the Mode i particle and having radius 
a, + aj lies inside the y'-th cylinder. We will denote that interior portion of the surface by Sy. 
Then if a search shows there are hy Mode j particles touching a particular Mode i particle 
identified by index number 4, then our best estimate of the total number of touching Mode j 
spheres, if the entire surface were accessible, would be 

rjij = h 0 4n( ai + ajf/Sij . VIII.C(l) (25n 

[00426] If the entire surface is accessible, then Sy = 47u(a i + aj) 2 and rjy = hy. Because of 
round-off error, a certain center-to-center error tolerance s is required to determine whether two 
spheres are indeed touching. The code considers them touching if their centers are within the 
distance (a, + aj)(l + e) from each other. Because the preferred implementation is in double 
precision, the default for £ is 10" 12 , but the user can specify any other desired e. Normally it 
should be a few orders of magnitude greater than machine precision. 

[00427] The spherical surface Sy lying inside a cylinder could be obtained from our 
earlier work on obtaining the spherical volume inside a cylinder by taking a derivative of the 
volume expressions with respect to the sphere radius. However, those expressions are 
complicated and have substantial bookkeeping associated with them. Because in this 
implementation there is an error tolerance, we may take a numerical derivative as follows: 

Sy = e x [V in (xih 9 yih Zfo, {a t + aj){\ + s); a c (j) 9 z botj , z top j) - 

V in (x ih , y ih , z ihi (at + aj)\ a c (j\ z both z top j)] . VIII.C(2) (252) 
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The 8 used in this expression should be kept very small even if the user chooses to be rather 
liberal in assigning a rather large error tolerance. One might, for example, assign this 8 to be 
either the user-defined (or default) error tolerance or 10" 6 , whichever is smaller. 



Finding Various Pack Statistics V'mdm^ Various Pack Statistics 

[00428] There are a number of vital statistics of the pack that will generally be of 
interest. If the user so specifies, then the following data is written to a user-named file. 

[00429] 1 . A recapitulation of the particle histogram. 

[00430] 2. A recapitulation of the values for control flags and parameters 
used in building the pack. 

[00431] 3. The total number of particles placed in the pack and the number 
placed from each mode. 

[00432] 4. The minimum and maximum radial and axial positions for the 
centers of all particles in the pack and the same for each mode. 

[00433] 5. Radial and axial particle density functions for each mode 
specified between the minimum and maximum positions specified above. 

[00434] 6. The number of particles placed after tunneling (if tunneling is 

permitted). 

[00435] 7. The number of particles hitting one of the catch nets. 

[00436] 8. The number of particles hitting one of the water levels. 

[00437] 9. The number of preplaced particles, if any. 
[00438] It is important to understanding in implementing the various aspects of the 
invention and in applying the related principles that the mathematical equations and expression 
presented herein may be assumed to have an infinite degree of precision and all real 
implementations of those equations result in answers with finite precision. Many programmers 
assume that the use of double precision, 6 4 bit 64-bit p recision, eliminates the need for close 
attention to this fact. Such an assumption with this problem typically will be incorrect. For 
example, in testing the implementation of several equations, cases were found where use of 
quadruple precision, 128 bit 128-bit p recision, was insufficient to obtain correct solutions. 
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These cases were avoided if in the solution process checks were made that allowed for tolerance 
bands. Every check for equality must include equality within a tolerance. The following 
tolerances are used within the preferred implementation. 

[00439] 1 . Comparison with z e ro. Comparison with zero. Numbers in the 
range -1.0E-20 < number < 1.0E-20 are treated as zero. 

[00440] 2. General compari s on for e quality. General comparison for 
equality. Any two non-zero numbers are within 1.0E-10 times the smaller of the two numbers 
are treated as equal. 

[00441] 3. General — angle — compari s on — for — equality. — General angle 
comparison for equality. Any two angles within 1.0E-6 radians are treated as equal. 

[00442] 4. Di s tance comparison. Distance comparison. Any two distances 
within 1.0E-6 times the smallest particle radius are treated as equal. 

[00443] 5. Distanc e s quared compari s on. Distance squared comparison. 
Any two numbers representing distance squared within 1 .OE-6 times the square of the smallest 
particle radius are treated as equal. 

[00444] Additional advantages and modifications will readily occur to those skilled in 
the art. Therefore, the invention in its broader aspects is not limited to the specific details, 
representative devices and methods, and illustrative examples shown and described. 
Accordingly, departures may be made from such details without departing from the spirit or 
scope of the general inventive concept as defined by the appended claims and their equivalents. 
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